AnĂ¡lisis de Series Temporales

IOT Analytics

library(dplyr)
library(tidyverse)
library(kableExtra)
library(DT) # DataTables



library(imputeTS) # Tratamiento de datos faltantes

library(plotly)   # GrĂ¡ficas
library(ggplot2)  # GrĂ¡ficas
library(ggfortify)
library(gridExtra)

library(forecast) # Tratamiento de series temporales
load(file="Datos/DF_Energia_GMinutos.RData")
load(file="Datos/DF_Energia_GHoras.RData")
load(file="Datos/DF_Energia_GDiaria.RData")
load(file="Datos/DF_Energia_GMensual.RData")

Objetivo

ProyecciĂ³n del consumo energĂ©tico a travĂ©s de series temporales.

Procedimiento

Estudiaremos los datos de consumo energético correspondientes a los años 2007 a 2009 a través de series temporales.

Haremos proyecciĂ³n de futuro para el año 2010 y compararemos los resultados con los datos reales de este Ăºltimo año (serie real).

Series que vamos a estudiar:

  • EvoluciĂ³n mensual del consumo para los años 2007-2009. PredicciĂ³n para el año 2010
  • EvoluciĂ³n diaria del consumo. Años 2007-2009. PredicciĂ³n para el año 2010
  • EvoluciĂ³n semanal del consumo. Años 2007-2009 . PredicciĂ³n para el año 2010

EvoluciĂ³n mensual del consumo. Años 2007-2009

Serie temporal y visualizaciĂ³n

Vamos a representar la serie por meses para cada submedidor y tambien para la energĂ­a global co granularidad mensual. En este caso, la frecuencia serĂ¡ 12, tenemos una observaciĂ³n por mes para los años 2007,8 y 2009.

seriesMensuales <- Granularidad_meses %>% filter(year != 2006 & year != 2010)
## Creamos un objeto TS
tsSM1_Mensual<- ts(seriesMensuales$Sub_metering_1, frequency=12, start=c(2007,1))
tsSM2_Mensual<- ts(seriesMensuales$Sub_metering_2, frequency=12, start=c(2007,1))
tsSM3_Mensual<- ts(seriesMensuales$Sub_metering_3, frequency=12, start=c(2007,1))
tsSM_GAP_Mensual<- ts(seriesMensuales$Global_active_power, frequency=12, start=c(2007,1))
## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Mensual, ts.colour = 'red', xlab = "Año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n mensual del consumo de energĂ­a de la cocina \n Años 2007-2009")# + scale_x_date(date_labels = '%b%')

autoplot(tsSM2_Mensual, xlab = "Año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n mensual del consumo de energĂ­a de la lavanderĂ­a \n Años 2007-2009")

autoplot(tsSM3_Mensual, ts.colour = 'blue', xlab = "Año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n mensual del consumo de energĂ­a del aire acondicionado y termo elĂ©ctrico \n Años 2007-2009")

autoplot(tsSM_GAP_Mensual, ts.colour = 'green', xlab = "Año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n mensual del consumo de energĂ­a \n Años 2007-2009")

Parece que los datos tienen estacionariedad (movimientos que se repiten cada año), esto ocurre para todas las variables.

PredicciĂ³n del año 2010 antes de descomponer la serie

Vamos a plicar regresiĂ³n lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrĂ¡tico medio.

SM1. Cocina

Modelo

fitSM1_Mensual <- tslm(tsSM1_Mensual ~ trend + season) 
summary(fitSM1_Mensual)

Call:
tslm(formula = tsSM1_Mensual ~ trend + season)

Residuals:
     Min       1Q   Median       3Q      Max 
-20522.7  -5190.9    215.7   6726.8  16702.8 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  66531.9     7294.5   9.121 4.21e-09 ***
trend         -173.1      200.7  -0.862 0.397393    
season2     -19137.6     9635.6  -1.986 0.059061 .  
season3      -1263.5     9641.8  -0.131 0.896879    
season4     -15856.4     9652.3  -1.643 0.114034    
season5      -6134.4     9666.9  -0.635 0.531966    
season6      -9988.0     9685.6  -1.031 0.313160    
season7     -27046.6     9708.4  -2.786 0.010505 *  
season8     -38674.8     9735.4  -3.973 0.000602 ***
season9     -10641.7     9766.3  -1.090 0.287159    
season10    -15040.7     9801.3  -1.535 0.138536    
season11     -6105.3     9840.3  -0.620 0.541071    
season12     -3378.2     9883.2  -0.342 0.735595    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11800 on 23 degrees of freedom
Multiple R-squared:  0.5796,    Adjusted R-squared:  0.3603 
F-statistic: 2.643 on 12 and 23 DF,  p-value: 0.02184

Resultados:

  • Valor del coeficiente de determinaciĂ³n: 0.5796. Es un valor bajo, el modelo Ăºnicamente explica el 58% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

PredicciĂ³n IC

PredicciĂ³n para los prĂ³ximos 20 dĂ­as

forecastfitSM1_Mensual <- forecast(fitSM1_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Mensual, ylab= "Watt-Hours", xlab="Fecha",
     main="PredicciĂ³n del consumo energĂ©tico del año 2010 \n  a travĂ©s de un modelo de regresiĂ³n")

SM2. Lavadero

Modelo

fitSM2_Mensual <- tslm(tsSM2_Mensual ~ trend + season) 
summary(fitSM2_Mensual)

Call:
tslm(formula = tsSM2_Mensual ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-26325  -4033   1223   6470  16873 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  83389.3     7617.6  10.947 1.35e-10 ***
trend         -898.9      209.6  -4.289 0.000274 ***
season2     -13475.5    10062.5  -1.339 0.193598    
season3       7508.4    10069.0   0.746 0.463406    
season4     -15023.1    10079.9  -1.490 0.149708    
season5      -7727.5    10095.2  -0.765 0.451777    
season6     -11683.0    10114.7  -1.155 0.259931    
season7     -19075.8    10138.6  -1.882 0.072621 .  
season8     -28408.6    10166.7  -2.794 0.010305 *  
season9      -9204.4    10199.1  -0.902 0.376160    
season10      4809.1    10235.6   0.470 0.642890    
season11     -1696.0    10276.3  -0.165 0.870356    
season12     -3191.5    10321.1  -0.309 0.759941    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 12320 on 23 degrees of freedom
Multiple R-squared:  0.6541,    Adjusted R-squared:  0.4736 
F-statistic: 3.624 on 12 and 23 DF,  p-value: 0.003891

Resultados:

  • Valor del coeficiente de determinaciĂ³n: 0.6541. Es un valor bajo, el modelo Ăºnicamente explica el 65% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

PredicciĂ³n IC

PredicciĂ³n para los prĂ³ximos 12 meses (año 2010)

forecastfitSM2_Mensual <- forecast(fitSM2_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Mensual, ylab= "Watt-Hours", xlab="Fecha",
     main="PredicciĂ³n del consumo energĂ©tico del año 2010 \n  a travĂ©s de un modelo de regresiĂ³n")

Se prevee una bajada del consumo

SM3. AC y termo

Modelo

fitSM3_Mensual <- tslm(tsSM3_Mensual ~ trend + season) 
summary(fitSM3_Mensual)

Call:
tslm(formula = tsSM3_Mensual ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-87112 -18443   4781  17647  82107 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   299382      23189  12.911 5.06e-12 ***
trend           1866        638   2.926  0.00761 ** 
season2       -50805      30631  -1.659  0.11077    
season3       -29034      30651  -0.947  0.35336    
season4       -64190      30684  -2.092  0.04768 *  
season5       -54165      30730  -1.763  0.09125 .  
season6       -85197      30790  -2.767  0.01097 *  
season7      -144850      30863  -4.693 1.00e-04 ***
season8      -171144      30948  -5.530 1.27e-05 ***
season9       -69624      31047  -2.243  0.03486 *  
season10      -52696      31158  -1.691  0.10430    
season11      -37281      31282  -1.192  0.24550    
season12        7092      31418   0.226  0.82341    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 37510 on 23 degrees of freedom
Multiple R-squared:  0.7577,    Adjusted R-squared:  0.6312 
F-statistic: 5.992 on 12 and 23 DF,  p-value: 0.0001238

Resultados:

  • Valor del coeficiente de determinaciĂ³n: 0.7577. Es un valor que podrĂ­a ser aceptable, el modelo Ăºnicamente explica el 76% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

PredicciĂ³n IC

PredicciĂ³n para los prĂ³ximos 12 meses (año 2010)

forecastfitSM3_Mensual <- forecast(fitSM3_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Mensual, ylab= "Watt-Hours", xlab="Fecha",
     main="PredicciĂ³n del consumo energĂ©tico del año 2010 \n  a travĂ©s de un modelo de regresiĂ³n")

Se prevee una bajada del consumo

ComparaciĂ³n de los coeficientes y resultados de cada modelo

DatosRealesSeriesMensuales <- Granularidad_meses %>% filter(year == 2010 )

RMSE_SM1_GranMensual<-accuracy(forecastfitSM1_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
R2_SM1_GranMensual<-0.5796
MASE_SM1_GranMensual<-accuracy(forecastfitSM1_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_1)[12]

RMSE_SM2_GranMensual<-accuracy(forecastfitSM2_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
R2_SM2_GranMensual<-0.6541
MASE_SM2_GranMensual<-accuracy(forecastfitSM2_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_2)[12]



RMSE_SM3_GranMensual<-accuracy(forecastfitSM3_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_3)[4]
R2_SM3_GranMensual<-0.7577
MASE_SM3_GranMensual<-accuracy(forecastfitSM3_Mensual  ,DatosRealesSeriesMensuales$Sub_metering_3)[12]

ResEvolMensual<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensual,RMSE_SM2_GranMensual,RMSE_SM3_GranMensual),
     # MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMensual),
      R2 = c(R2_SM1_GranMensual,R2_SM2_GranMensual,R2_SM3_GranMensual)
)

ResEvolMensual %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE R2
Cocina 11186.74 0.5796
Lavadero 10939.60 0.6541
AC y TermoE 57300.05 0.7577
# MASE: Mean absolute error
# RMSE: root mean scuare error. Diferencia entre los valores predichos y los observados.

El RMSE del modelo con mejor coeficiente de determinaciĂ³n es el mĂ¡s alto. Pero no hay que olvidar que el consumo elĂ©ctrico para el submedidor correspondiente al AC y termo elĂ©ctrico es el que mĂ¡s energĂ­a consume de todos con una diferencia muy grande.

Forecasting descomponiendo la serie (asĂ­ ok)

Descomponer la serie temporal: identificar y calcular las componentes tendencia, estacionalidad y ruido.

  • Tendencia de una serie temporal: Consiste en la evoluciĂ³n a largo plazo de la serie. En nuestro caso, puede haber una tendencia al aumento del consumo energĂ©tico o una tendencia decreciente.
  • Estacionalidad: Movimientos oscilatorios alrededor de la tendencia que se repiten de forma periĂ³dica y con amplitud inferior al año
  • Componente irregular: es aleatoria y no se puede predecir (ruido)

DescomposiciĂ³n de la serie y visualizaciĂ³n

Cocina: Submetering 1

Componentes_SM1_MensualSTL<- tsSM1_Mensual %>%stl( t.window=12,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_MensualSTL, ts.colour = 'red', main="DescomposiciĂ³n de componentes")

Se aprecia una gran influencia de la componente estacional: oscilaciones que se repiten cada año, por ejemplo, en el mes de agosto, el consumo energético de la cocina decrece consuderablemente y es mucho menor que el resto de meses. Por el contrario, es en el mes de Junio cuando se registra el mayor consumo energético de los aparatos submeterizados de la cocina.

Tambien se aprecia una tendencia creciente de consumo, cada año se consume mĂ¡s.

Remainder no es un muelle, grĂ¡ficamente, parece que la componente irregular no es del todo aleatoria.

plot.ts(tsSM1_Mensual, col = 'gray',
         xlab = "Mes del año", 
        ylab = "EnergĂ­a (Varios-hora)", 
        main = "EvoluciĂ³n mensual del consumo energĂ©tico de la cocina \n y la tendencia del consumo")
lines(Componentes_SM1_MensualSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Como ya habíamos comentado, existe una tendencia creciente del consumo energético de la cocina.

Lavadero: Submetering 2

Componentes_SM2_MensualSTL<- tsSM2_Mensual %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_MensualSTL, ts.colour = 'blue')

Por el contrario que en la cocina, se observa como existe una considerable tendencia decreciente del consumo en el lavadero.

Tambien se observa claramente la estacionalidad, donde podemos apreciar que en los meses de verano el consumo energético es mucho menor que el resto de meses del año.

Parece que la componente irregular podrĂ­a ser aleatoria, sin embargo, esto es Ăºnicamente a modo grĂ¡fico y por tanto deberĂ­amos comprobarlo con un test de hipĂ³tesis.

plot.ts(tsSM2_Mensual, col = 'gray', main="EvoluciĂ³n mensual del consumo energĂ©tico del lavadero")
lines(Componentes_SM2_MensualSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2
      )

Termo eléctrico y aire acondicionado: Submetering 3

tsSM3_Mensual %>%
  stl( t.window=12,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", t.window = 12, robust = TRUE)

 Time.series components:
    seasonal              trend            remainder        
 Min.   :-107286.17   Min.   :236737.3   Min.   :-98598.49  
 1st Qu.: -12301.07   1st Qu.:257514.6   1st Qu.:-14547.16  
 Median :   8318.49   Median :277783.5   Median :  -493.40  
 Mean   :      0.00   Mean   :274647.8   Mean   : -3392.79  
 3rd Qu.:  27345.56   3rd Qu.:287394.0   3rd Qu.:  9934.48  
 Max.   :  72770.26   Max.   :312023.9   Max.   : 91583.87  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     39647        29879     24482         82030
   %  48.3         36.4      29.8         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8138  0.9452  0.8280  0.9944  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 361 12 13
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 37 2 2
 $ inner: int 1
 $ outer: int 15
Componentes_SM3_MensualSTL<- tsSM3_Mensual %>%
  stl( t.window=12, s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_MensualSTL, ts.colour = 'blue')

Tendencia creciente al consumo energĂ©tico. Se observa fuerte componente estacional, donde el consumo energĂ©tico baja considerablemente en los meses de verano, y alcanza niveles mĂ¡ximos en los Ăºltimos meses del año (Noviembre y Diciembre)

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente

plot.ts(tsSM3_Mensual, col = 'gray',
        main="EvoluciĂ³n del consumo energĂ©tico mensual del aire acondicionado y termo elĂ©ctrico")
lines(Componentes_SM3_MensualSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

PredicciĂ³n con Holt-Winters (suavizado exponencial) para el año 2010

Vamos a aplicar suavizado exponencial para realizar un pronĂ³stico del comportamiento de consumo energĂ©tico. Veremos esta predicciĂ³n para datos:

  • Sin estacionalidad
  • Con estacionalidad, para ver si se mantiene o no la tendencia.

Tendencia de predicciĂ³n (crecimiento o decrecimiento del consumo) y predicciĂ³n teniendo en cuenta las componentes

EvoluciĂ³n mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina
tsSM1_Mensual_detrended <- tsSM1_Mensual - Componentes_SM1_MensualSTL$time.series[,1]
# Elimino la componente estacional
plot.ts(tsSM1_Mensual_detrended, 
        ylab = "Energía (Vatios-hora)", xlab="Mes del año", #ylim=c(0,3000),
        main = 'EvoluciĂ³n del consumo energĂ©tico de la cocina sin estacionalidad',
        cex.main = 0.85, col="gray")

Vamos a estudiar como se comporta la componente irregular.

par(mfrow = c(1,2))
plot(Componentes_SM1_MensualSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_MensualSTL$time.series[,3])
qqline(Componentes_SM1_MensualSTL$time.series[,3])

shapiro.test(Componentes_SM1_MensualSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM1_MensualSTL$time.series[, 3]
W = 0.75719, p-value = 2.554e-06

Los resĂ­duos no son solo ruido. AĂºn hay estacionalidad.

Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_Mensual_detrended ,s.window = "periodic",robust = TRUE),colour = 'red', main="DescomposiciĂ³n de la serie desestacionalizada")

PodrĂ­amos a sumir que se ha eliminado la estacionalidad, ya que los resĂ­duos se comportan de manera irregular.

Suavizado exponencial de HoltWinters
  • Primera iteracciĂ³n: Solo con \(\alpha\), voy a averiguar en que margen se va a mover la predicciĂ³n
tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , seasonal = "additive",
#   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM1_HW_Mensual$fitted)
      xhat           level           trend            season     
 Min.   : 7683   Min.   :13713   Min.   :-376.9   Min.   :-7702  
 1st Qu.:45514   1st Qu.:50966   1st Qu.:-376.9   1st Qu.:-5695  
 Median :49998   Median :52578   Median :-376.9   Median :-1479  
 Mean   :49432   Mean   :49809   Mean   :-376.9   Mean   :    0  
 3rd Qu.:57209   3rd Qu.:58150   3rd Qu.:-376.9   3rd Qu.: 1006  
 Max.   :88889   Max.   :66538   Max.   :-376.9   Max.   :25965  
# optimiza
plot(tsSM1_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronĂ³stico de consumo serĂ¡ en el rango: 7683-88889 vatios-hora de energĂ­a consumida.

  • Segunda iteracciĂ³n: Añadiendo \(\beta\) podrĂ© averiguar la tendencia del consumo para las predicciones
tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM1_HW_Mensual$fitted)
      xhat            level           trend              season     
 Min.   :-19652   Min.   :13713   Min.   :-36810.3   Min.   :-7702  
 1st Qu.: 44754   1st Qu.:50966   1st Qu.: -3854.5   1st Qu.:-5695  
 Median : 50051   Median :52578   Median :  -188.6   Median :-1479  
 Mean   : 50152   Mean   :49809   Mean   :   343.0   Mean   :    0  
 3rd Qu.: 69761   3rd Qu.:58150   3rd Qu.:  5871.7   3rd Qu.: 1006  
 Max.   :101263   Max.   :66538   Max.   : 42426.8   Max.   :25965  
# optimiza
plot(tsSM1_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Tendencia al consumo decreciente, es decir, se pronosticarĂ¡ un descenso del consumo par meses siguientes.

  • Tercera iteracciĂ³n: Si añado \(\gamma\) podrĂ© ver la estacionalidad.
tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
 # gamma=TRUE,
   alpha=TRUE
)

summary(tsSM1_HW_Mensual$fitted)
      xhat            level           trend              season     
 Min.   :-19652   Min.   :13713   Min.   :-36810.3   Min.   :-7702  
 1st Qu.: 44754   1st Qu.:50966   1st Qu.: -3854.5   1st Qu.:-5695  
 Median : 50051   Median :52578   Median :  -188.6   Median :-1479  
 Mean   : 50152   Mean   :49809   Mean   :   343.0   Mean   :    0  
 3rd Qu.: 69761   3rd Qu.:58150   3rd Qu.:  5871.7   3rd Qu.: 1006  
 Max.   :101263   Max.   :66538   Max.   : 42426.8   Max.   :25965  
# optimiza
plot(tsSM1_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

# Rstudio optimiza los coeficientes
tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , seasonal = "additive",
 # beta=TRUE, 
 # gamma=TRUE,
 #  alpha=TRUE
)

summary(tsSM1_HW_Mensual$fitted)
      xhat           level           trend            season       
 Min.   :41445   Min.   :46024   Min.   :-376.9   Min.   :-7930.2  
 1st Qu.:44274   1st Qu.:48191   1st Qu.:-376.9   1st Qu.:-5667.3  
 Median :47013   Median :50358   Median :-376.9   Median :-1903.9  
 Mean   :49704   Mean   :50358   Mean   :-376.9   Mean   : -276.8  
 3rd Qu.:51922   3rd Qu.:52525   3rd Qu.:-376.9   3rd Qu.: 1093.8  
 Max.   :76135   Max.   :54692   Max.   :-376.9   Max.   :25965.5  
# optimiza
plot(tsSM1_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_for <- forecast(tsSM1_HW_Mensual, h=12)
plot(tsSM1_HW_Mensual_for ,
     ylab= "Vatios-Hora",
     xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="PredicciĂ³n del consumo energĂ©tico de la cocina \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forC <- forecast(tsSM1_HW_Mensual, h=12
                                  , 
                                  level=c(10,25)
                                  )
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha",
     main="PronĂ³stico del consumo de la cocina. Año 2010 \n (sin estacionalidad)", start(2010))

Lavadero
tsSM2_Mensual_detrended <- tsSM2_Mensual - Componentes_SM2_MensualSTL$time.series[,2]
plot.ts(tsSM2_Mensual_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", col="gray",
        main = 'Variabilidad del consumo energético del lavadero \n para datos sin estacionalidad',
        cex.main = 0.85)

Vamos a estudiar como se comporta la componente irregular.

par(mfrow = c(1,2))
plot(Componentes_SM2_MensualSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_MensualSTL$time.series[,3])
qqline(Componentes_SM2_MensualSTL$time.series[,3])

shapiro.test(Componentes_SM2_MensualSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM2_MensualSTL$time.series[, 3]
W = 0.85783, p-value = 0.0002898

Los resĂ­duos no siguen una distribuciĂ³n normal. La descomposiciĂ³n no es buena, por tanto, tenemos que vovler a descomponer la serie.

## Volvemos a descomponer la serie
plot(stl(tsSM2_Mensual_detrended ,s.window = "periodic",
         robust = TRUE,
         t.window = 12))

Ahora si que parece que los resíduos se comportan de manera aleatoria. Comportamiento de consumo: Hay una tendencia decreciente al consumo los 5 primeros meses del año 2007, sin embargo, posteriormente este comportamiento cambia y existe una tendencia al consumo hasta principios del año 2008, donde vuelve a cambiar la tendencia y pasa a haber un comportamiento decreciente hasya finales del año 2008. Al parecer, desde final del año 2008 se ha observado una tendencia creciente del consumo energético del lavadero.

Suavizado exponencial de HoltWinters
  • Primera iteracciĂ³n: Solo con \(\alpha\), voy a averiguar en que margen se va a mover la predicciĂ³n
tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , seasonal = "additive",
#   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM2_HW_Mensual$fitted)
      xhat              level            trend           season      
 Min.   :-37423.2   Min.   :-21787   Min.   :237.4   Min.   :-22162  
 1st Qu.: -9524.1   1st Qu.: -7651   1st Qu.:237.4   1st Qu.: -8854  
 Median : -1669.0   Median :  1753   Median :237.4   Median :  3244  
 Mean   :  -131.6   Mean   :  -369   Mean   :237.4   Mean   :     0  
 3rd Qu.: 11559.8   3rd Qu.:  3709   3rd Qu.:237.4   3rd Qu.:  8595  
 Max.   : 27452.5   Max.   : 32063   Max.   :237.4   Max.   : 17041  
# optimiza
plot(tsSM2_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronĂ³stico de consumo serĂ¡ en el rango: 0-27452.5 vatios-hora de energĂ­a consumida, ya que no puede haber un consumo energĂ©tico negativo.

  • Segunda iteracciĂ³n: Añadiendo \(\beta\) podrĂ© averiguar la tendencia del consumo para las predicciones
tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM2_HW_Mensual$fitted)
      xhat              level            trend              season      
 Min.   :-62015.7   Min.   :-21787   Min.   :-27388.9   Min.   :-22162  
 1st Qu.:-13466.1   1st Qu.: -7651   1st Qu.: -6777.2   1st Qu.: -8854  
 Median :  2156.0   Median :  1753   Median :  -158.7   Median :  3244  
 Mean   :  -144.9   Mean   :  -369   Mean   :   224.2   Mean   :     0  
 3rd Qu.: 16152.5   3rd Qu.:  3709   3rd Qu.:  7242.4   3rd Qu.:  8595  
 Max.   : 73833.0   Max.   : 32063   Max.   : 53850.3   Max.   : 17041  
# optimiza
plot(tsSM2_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Tendencia al consumo cambiante como ya hemos comentado antes

  • Tercera iteracciĂ³n: Si añado \(\gamma\) podrĂ© ver la estacionalidad.
tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)

summary(tsSM2_HW_Mensual$fitted)
      xhat              level            trend              season      
 Min.   :-62015.7   Min.   :-21787   Min.   :-27388.9   Min.   :-22162  
 1st Qu.:-13466.1   1st Qu.: -7651   1st Qu.: -6777.2   1st Qu.: -8854  
 Median :  2156.0   Median :  1753   Median :  -158.7   Median :  3244  
 Mean   :  -144.9   Mean   :  -369   Mean   :   224.2   Mean   :     0  
 3rd Qu.: 16152.5   3rd Qu.:  3709   3rd Qu.:  7242.4   3rd Qu.:  8595  
 Max.   : 73833.0   Max.   : 32063   Max.   : 53850.3   Max.   : 17041  
# optimiza
plot(tsSM2_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

# Rstudio optimiza los coeficientes
tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE,
 #  alpha=TRUE
)

summary(tsSM2_HW_Mensual$fitted)
      xhat               level             trend           season        
 Min.   :-28433.55   Min.   :-1983.2   Min.   :234.8   Min.   :-31197.1  
 1st Qu.: -7760.62   1st Qu.: -612.7   1st Qu.:236.3   1st Qu.: -9057.5  
 Median :  2416.46   Median :  755.1   Median :237.5   Median :  1068.6  
 Mean   :    86.66   Mean   :  751.1   Mean   :237.2   Mean   :  -901.7  
 3rd Qu.:  9550.00   3rd Qu.: 2117.5   3rd Qu.:238.2   3rd Qu.:  8271.5  
 Max.   : 18504.81   Max.   : 3472.5   Max.   :239.1   Max.   : 17041.2  
# optimiza
plot(tsSM2_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_for <- forecast(tsSM2_HW_Mensual, h=12)
plot(tsSM2_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n del consumo energĂ©tico del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forC <- forecast(tsSM2_HW_Mensual, h=12, level=95)
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n del consumo energĂ©tico del lavadero. Año 2010", start(2010))

AC y termo eléctrico
tsSM3_Mensual_detrended <- tsSM3_Mensual - Componentes_SM3_MensualSTL$time.series[,2]
plot.ts(tsSM3_Mensual_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'EvoluciĂ³n el consumo energĂ©tico del aire acondicionado y termo ele´ctrico \n Datos sin estacionalidad',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM3_MensualSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_MensualSTL$time.series[,3])
qqline(Componentes_SM3_MensualSTL$time.series[,3])

shapiro.test(Componentes_SM3_MensualSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM3_MensualSTL$time.series[, 3]
W = 0.87556, p-value = 0.0007821

Los resĂ­duos no siguen una distribuciĂ³n normal. La descomposiciĂ³n no es buena.

Hay que volver a descomponer la serie.

## Volvemos a descomponer la serie
plot(stl(tsSM3_Mensual_detrended ,s.window = "periodic",robust = TRUE))

Parece que la componente irregular es aleatoria. La tendencia del consumo energético es creciente, aunque durante el año 2008 hubo períodos con tendencia decreciente, al igual que los 5 primeros meses del año 2007.

Suavizado exponencial de HoltWinters
  • Primera iteracciĂ³n: Solo con \(\alpha\), voy a averiguar en que margen se va a mover la predicciĂ³n
tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended , seasonal = "additive",
#   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM3_HW_Mensual$fitted)
      xhat             level             trend          season      
 Min.   :-209442   Min.   :-181400   Min.   :-635   Min.   :-92079  
 1st Qu.:  -9333   1st Qu.:  -8991   1st Qu.:-635   1st Qu.:-23534  
 Median :   9306   Median :   4859   Median :-635   Median :   699  
 Mean   :  -5753   Mean   :  -5118   Mean   :-635   Mean   :     0  
 3rd Qu.:  31822   3rd Qu.:  20231   3rd Qu.:-635   3rd Qu.: 17164  
 Max.   :  80208   Max.   :  45472   Max.   :-635   Max.   : 86727  
# optimiza
plot(tsSM3_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronĂ³stico de consumo serĂ¡ en el rango: 0-80208 vatios-hora de energĂ­a consumida, ya que no puede haber un consumo energĂ©tico negativo.

  • Segunda iteracciĂ³n: Añadiendo \(\beta\) podrĂ© averiguar la tendencia del consumo para las predicciones
tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
#   gamma=TRUE,
   alpha=TRUE
)

summary(tsSM3_HW_Mensual$fitted)
      xhat             level             trend               season      
 Min.   :-422399   Min.   :-181400   Min.   :-213592.3   Min.   :-92079  
 1st Qu.: -27389   1st Qu.:  -8991   1st Qu.: -26879.7   1st Qu.:-23534  
 Median :  11408   Median :   4859   Median :  -1375.6   Median :   699  
 Mean   :  -5558   Mean   :  -5118   Mean   :   -439.2   Mean   :     0  
 3rd Qu.:  53303   3rd Qu.:  20231   3rd Qu.:  15913.9   3rd Qu.: 17164  
 Max.   : 245407   Max.   :  45472   Max.   : 213552.8   Max.   : 86727  
# optimiza
plot(tsSM3_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Tendencia al consumo cambiante como ya hemos comentado antes.

  • Tercera iteracciĂ³n: Si añado \(\gamma\) podrĂ© ver la estacionalidad.
tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended , seasonal = "additive",
   beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)

summary(tsSM3_HW_Mensual$fitted)
      xhat             level             trend               season      
 Min.   :-422399   Min.   :-181400   Min.   :-213592.3   Min.   :-92079  
 1st Qu.: -27389   1st Qu.:  -8991   1st Qu.: -26879.7   1st Qu.:-23534  
 Median :  11408   Median :   4859   Median :  -1375.6   Median :   699  
 Mean   :  -5558   Mean   :  -5118   Mean   :   -439.2   Mean   :     0  
 3rd Qu.:  53303   3rd Qu.:  20231   3rd Qu.:  15913.9   3rd Qu.: 17164  
 Max.   : 245407   Max.   :  45472   Max.   : 213552.8   Max.   : 86727  
# optimiza
plot(tsSM3_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El ajuste podrĂ­a ser mejorable, ya que hay gran influencia de los valores extremos.

## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended ,
                              # tsSM3_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Mensual,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_for <- forecast(tsSM3_HW_Mensual, h=12)
plot(tsSM3_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n del consumo energĂ©tico del AC y termo elĂ©ctrico \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forC <- forecast(tsSM3_HW_Mensual, h=12, level=95)
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n el consumo energĂ©tico del lavadero y termo elĂ©ctrico para el año 2010", start(2010),
     ylim=c(0,80000))

ComparaciĂ³n de los coeficientes y resultados de cada predicciĂ³n. Datos sin estacionalidad

En las grĂ¡ficas de predicciĂ³n, hemos observado que el ajuste no mantiene la tendencia de consumo energĂ©tico de la serie.

RMSE_SM1_GranMensualFor<-accuracy(tsSM1_HW_Mensual_forC   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualFor<-accuracy(tsSM2_HW_Mensual_forC   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualFor<-accuracy(tsSM3_HW_Mensual_forC   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2_SinEstacionalidadConTendencia <-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualFor,
               RMSE_SM2_GranMensualFor,
               RMSE_SM3_GranMensualFor)
)
ResEvolMensual2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 14278.84
Lavadero 44294.27
AC y TermoE 325349.90

EvoluciĂ³n mensual del consumo (con estacionalidad). Años 2007-2009

Vamos a ver una predicciĂ³n del consumo enegĂ©tico de cada submedidor con el mĂ©todo de HW. En este caso, no eliminaremos la estacionalidad, para ver asĂ­ si se mantiene la tendencia de consumo.

Cocina
tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)

## Volvemos a descomponer la serie
plot(stl(tsSM1_Mensual_AdjustedII,s.window = "periodic",robust = TRUE,t.window = 12))

Tendencia creciente del consumo.

Suavizado exponencial de HoltWinters
  • Primera interacciĂ³n: para observar el rango de la predicciĂ³n.
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, alpha=TRUE,
                                #  beta=TRUE, gamma=TRUE
                                  )
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

- Segunda interacciĂ³n: añado \(\beta\) para observar la tendencia del consumo:

## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, alpha=TRUE,
                                  beta=TRUE # , gamma=TRUE
                                  )
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

- Tercera iteracciĂ³n: Añado \(\gamma\)

## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII
                                  , 
                                  alpha=TRUE,
                                  beta=TRUE  , gamma=TRUE
                                  )
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, 
                                 #  alpha=TRUE,
                                   beta=TRUE  ,
                                  gamma=TRUE
                                  )
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Si dejamos que la propia funciĂ³n optimice los valores de beta y gamma, la predicciĂ³n se ajusta mejor y por tanto, obtendremos mejores resultados.

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="PredicciĂ³n del consumo energĂ©tico de la cocina \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=12, level=95)
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha", main="PronĂ³stico para el año 2010 (con estacionalidad)", start(2010))

Lavadero
tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)

## Volvemos a descomponer la serie
plot(stl(tsSM2_Mensual_AdjustedII,s.window = "periodic",
         robust = TRUE,t.window = 12))

LA tendencia de consumo es claramente creciente, y parece que la componente irregular se ha eliminado.

Suavizado exponencial de HoltWinters
  • Primera interacciĂ³n: componente alpha para ver el rango de predicciones
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, 
                                  # beta=TRUE, 
                                  # gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM2_HW_MensualII$fitted)
      xhat           level           trend           season      
 Min.   :10670   Min.   :27252   Min.   :-1121   Min.   :-22206  
 1st Qu.:41144   1st Qu.:44247   1st Qu.:-1121   1st Qu.: -8213  
 Median :56659   Median :56542   Median :-1121   Median :  3157  
 Mean   :52598   Mean   :53719   Mean   :-1121   Mean   :     0  
 3rd Qu.:64844   3rd Qu.:62657   3rd Qu.:-1121   3rd Qu.:  8459  
 Max.   :80481   Max.   :80053   Max.   :-1121   Max.   : 17218  
plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Parece que la predicciĂ³n del consumo energĂ©tico del lavadero estarĂ¡ entre 10670 y 80481 vatios-hora de energĂ­a.

  • Segunda interacciĂ³m: añado beta para ver la tendencia del consumo:
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, 
                                   beta=TRUE, 
                                  # gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM2_HW_MensualII$fitted)
      xhat            level           trend              season      
 Min.   :-12375   Min.   :27252   Min.   :-27249.7   Min.   :-22206  
 1st Qu.: 39497   1st Qu.:44247   1st Qu.: -7920.3   1st Qu.: -8213  
 Median : 56536   Median :56542   Median : -2034.5   Median :  3157  
 Mean   : 52768   Mean   :53719   Mean   :  -951.5   Mean   :     0  
 3rd Qu.: 66700   3rd Qu.:62657   3rd Qu.:  6611.9   3rd Qu.:  8459  
 Max.   :121207   Max.   :80053   Max.   : 52800.2   Max.   : 17218  
plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Como ya habĂ­amos avanzado, tendencia creciente

  • Tercera interacciĂ³n: añadimos gamma
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, 
                                   beta=TRUE, 
                                   gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM2_HW_MensualII$fitted)
      xhat            level           trend              season      
 Min.   :-12375   Min.   :27252   Min.   :-27249.7   Min.   :-22206  
 1st Qu.: 39497   1st Qu.:44247   1st Qu.: -7920.3   1st Qu.: -8213  
 Median : 56536   Median :56542   Median : -2034.5   Median :  3157  
 Mean   : 52768   Mean   :53719   Mean   :  -951.5   Mean   :     0  
 3rd Qu.: 66700   3rd Qu.:62657   3rd Qu.:  6611.9   3rd Qu.:  8459  
 Max.   :121207   Max.   :80053   Max.   : 52800.2   Max.   : 17218  
plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, 
                                   beta=TRUE, 
                                   gamma=TRUE #,
                                 #  alpha = TRUE
                                  )
PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=12)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n del consumo energĂ©tico del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=12, level=95)
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n para el año 2010", start(2010))

AC y termo eléctrico
tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)

## Volvemos a descomponer la serie
plot(stl(tsSM3_Mensual_AdjustedII,s.window = "periodic",robust = TRUE, 
         t.window = 12))

Tendencia creciente del consumo de este submedidor.

Suavizado exponencial de HoltWinters
  • Primera interacciĂ³n: alpha para ver el margen de predicciĂ³n
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII,
                                #  beta=TRUE, 
                                #  gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM3_HW_MensualII$fitted)
      xhat            level            trend          season       
 Min.   : 80833   Min.   :111338   Min.   :2186   Min.   :-104479  
 1st Qu.:262195   1st Qu.:268973   1st Qu.:2186   1st Qu.: -24576  
 Median :295777   Median :280868   Median :2186   Median :   6764  
 Mean   :281117   Mean   :278931   Mean   :2186   Mean   :      0  
 3rd Qu.:319558   3rd Qu.:307632   3rd Qu.:2186   3rd Qu.:  21778  
 Max.   :392649   Max.   :331536   Max.   :2186   Max.   :  86499  
plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Rrango de predicciĂ³nes: 80833-392694 vatios hora.

  • Segunda interacciĂ³n: añado beta para ver la tendencia del consumo
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII,
                                  beta=TRUE, 
                                #  gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM3_HW_MensualII$fitted)
      xhat             level            trend             season       
 Min.   :-141552   Min.   :111338   Min.   :-220198   Min.   :-104479  
 1st Qu.: 268345   1st Qu.:268973   1st Qu.: -28339   1st Qu.: -24576  
 Median : 293751   Median :280868   Median :   1066   Median :   6764  
 Mean   : 281347   Mean   :278931   Mean   :   2417   Mean   :      0  
 3rd Qu.: 352536   3rd Qu.:307632   3rd Qu.:  23657   3rd Qu.:  21778  
 Max.   : 524064   Max.   :331536   Max.   : 207838   Max.   :  86499  
plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Parece que los valores extremos influirĂ¡n en las predicciones

  • Tercera interacciĂ³n
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII,
                                  beta=TRUE, 
                                  gamma=TRUE,
                                  alpha = TRUE
                                  )
summary(tsSM3_HW_MensualII$fitted)
      xhat             level            trend             season       
 Min.   :-141552   Min.   :111338   Min.   :-220198   Min.   :-104479  
 1st Qu.: 268345   1st Qu.:268973   1st Qu.: -28339   1st Qu.: -24576  
 Median : 293751   Median :280868   Median :   1066   Median :   6764  
 Mean   : 281347   Mean   :278931   Mean   :   2417   Mean   :      0  
 3rd Qu.: 352536   3rd Qu.:307632   3rd Qu.:  23657   3rd Qu.:  21778  
 Max.   : 524064   Max.   :331536   Max.   : 207838   Max.   :  86499  
plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=12)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n del consumo energĂ©tico del AC y termo elĂ©ctrico \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=12, level=95)
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n del consumo energĂ©tico del aire acondicionado \n  y termo elĂ©ctrico  para el año 2010", start(2010), 
     ylim=c(0,400000))

ComparaciĂ³n de los coeficientes y resultados de cada predicciĂ³n. Datos con estacionalidad
RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualForII,
               RMSE_SM2_GranMensualForII,
               RMSE_SM3_GranMensualForII)
   
    
)

ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 16362.16
Lavadero 27223.90
AC y TermoE 83037.90

ComparaciĂ³n de resultados. PredicciĂ³n antes y despues de eliminar estacionalidad

Comparacion<-cbind.data.frame(
  ResEvolMensual2_SinEstacionalidadConTendencia,
  ResEvolMensual2SinTendenciaConEstacionalidad[,2])

colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE Con Estacionalidad RMSE Sin estacionalidad
Cocina 14278.84 16362.16
Lavadero 44294.27 27223.90
AC y TermoE 325349.90 83037.90

Resultados mĂ¡s fiables y con menor error trĂ¡s haber eliminado la estacionalidad de la serie. TambiĂ©n asĂ­ podemos ver la tendencia.

ComparaciĂ³n de resultados. PredicciĂ³n antes y despues de eliminar la tendencia (descomponer)

En las grĂ¡ficas de predicciĂ³n de suaviazado exponencial, hemos observado que el ajuste cuando eliminamos la estacionalidad, no mantiene la tendencia de consumo energĂ©tico de la serie. Sin embargo, al mantener la estacionalidad se puede apreciar como las predicciones ajustan la tendencia de la serie original

# PredicciĂ³n antes y despues de eliminar estacionalidad:
Comparacion
   Submedidor RMSE Con Estacionalidad RMSE Sin estacionalidad
1      Cocina                14278.84                16362.16
2    Lavadero                44294.27                27223.90
3 AC y TermoE               325349.90                83037.90
# Sin descomponer


ComparacionII<-cbind.data.frame(ResEvolMensual,
                              Comparacion[-1]
                              )
colnames(ComparacionII) <- c("Submedidor","RMSE","R^2","Con estacionalidad","Sin estacionalidad")
ComparacionII %>% kable(booktabs=TRUE) %>% 
  kable_styling(latex_options = "striped") %>% 
  add_header_above( c("", "Serie sin descomponer" = 2, "Descomponiendo la serie" = 2)) %>% 
  add_header_above( c("", "RMSE" = 4))
RMSE
Serie sin descomponer
Descomponiendo la serie
Submedidor RMSE R^2 Con estacionalidad Sin estacionalidad
Cocina 11186.74 0.5796 14278.84 16362.16
Lavadero 10939.60 0.6541 44294.27 27223.90
AC y TermoE 57300.05 0.7577 325349.90 83037.90

Por tanto, nos quedaremos con los modelos correspondientes a la serie sin estacionalidad, que son los que mejores resultados nos ha dado.

PronĂ³stico de consumo para los prĂ³ximos 12 meses

Cocina

# DefiniciĂ³n de serie temporal
seriesMensualesCompletas <- Granularidad_meses 
## Creamos un objeto TS
tsSM1_Mensual<- ts(seriesMensualesCompletas$Sub_metering_1, frequency=12, start=c(2006,12))
tsSM2_Mensual<- ts(seriesMensualesCompletas$Sub_metering_2, frequency=12, start=c(2006,12))
tsSM3_Mensual<- ts(seriesMensualesCompletas$Sub_metering_3, frequency=12, start=c(2006,12))
tsSM_GAP_Mensual<- ts(seriesMensualesCompletas$Global_active_power, frequency=12, start=c(2006,12))
# DescomposiciĂ³n de la serie
Componentes_SM1_MensualSTL<-  tsSM1_Mensual %>%
  stl( t.window=12,s.window="periodic",  robust=TRUE)
# EliminaciĂ³n de la estacionalidad
tsSM1_Mensual_detrended <- tsSM1_Mensual - Componentes_SM1_MensualSTL$time.series[,1]
# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended , 
                                seasonal = "additive",
 # beta=TRUE, 
 # gamma=TRUE,
 #  alpha=TRUE
)
# PronĂ³stico para los siguientes 13 meses 
PrediccionFinalCocinaAnual <- forecast(tsSM1_HW_Mensual, h=13)
## Plot only the forecasted area
plot(PrediccionFinalCocinaAnual, ylab= "Vatios-Hora", xlab="Fecha", main="PredicciĂ³n el consumo energĂ©tico \n  de la cocina para los prĂ³ximos 12 meses", start(2010,12),
     ylim=c(0,80000))

Las predicciones son las siguientes:

PrediccionFinalCocinaAnual[["mean"]] 
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2010                                                                        
2011 41945.83 39455.73 37523.93 44392.00 45300.17 34904.71 39853.79 39158.80
          Sep      Oct      Nov      Dec
2010                            65059.25
2011 42425.45 39842.72 39116.19 62084.64

Lavadero

# DescomposiciĂ³n de la serie
Componentes_SM2_MensualSTL<-  tsSM2_Mensual %>%
  stl( t.window=12,s.window="periodic",  robust=TRUE)
# EliminaciĂ³n de la estacionalidad
tsSM2_Mensual_detrended <- tsSM2_Mensual - Componentes_SM2_MensualSTL$time.series[,1]
# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended , 
                                seasonal = "additive",
 # beta=TRUE, 
 # gamma=TRUE,
 #  alpha=TRUE
)
# PronĂ³stico para los siguientes 13 meses 
PrediccionFinalLavaderoAnual <- forecast(tsSM2_HW_Mensual, h=13)
## Plot only the forecasted area
plot(PrediccionFinalLavaderoAnual, ylab= "Vatios-Hora", xlab="Fecha", main="PredicciĂ³n el consumo energĂ©tico \n  del lavadero para los prĂ³ximos 12 meses", start(2010,12),
     ylim=c(0,80000))

Las predicciones son las siguientes:

PrediccionFinalLavaderoAnual[["mean"]] 
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2010                                                                        
2011 34865.05 38962.01 33876.42 30611.39 37750.86 30764.14 35867.64 33948.36
          Sep      Oct      Nov      Dec
2010                            56335.97
2011 35359.50 31508.32 31995.30 45813.51

Aire acondicionado y termo eléctrico

# DescomposiciĂ³n de la serie
Componentes_SM3_MensualSTL<-  tsSM3_Mensual %>%
  stl( t.window=12,s.window="periodic",  robust=TRUE)
# EliminaciĂ³n de la estacionalidad
tsSM3_Mensual_detrended <- tsSM3_Mensual - Componentes_SM3_MensualSTL$time.series[,1]
# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended , 
                                seasonal = "additive",
 # beta=TRUE, 
 # gamma=TRUE,
 #  alpha=TRUE
)
# PronĂ³stico para los siguientes 13 meses 
PrediccionFinalAcTAnual <- forecast(tsSM3_HW_Mensual, h=13)
## Plot only the forecasted area
plot(PrediccionFinalAcTAnual, ylab= "Vatios-Hora", xlab="Fecha", main="PredicciĂ³n el consumo energĂ©tico \n  del aire y termo elĂ©ctricopara los prĂ³ximos 12 meses", start(2010,12))

Las predicciones son las siguientes:

PrediccionFinalLavaderoAnual[["mean"]] 
          Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
2010                                                                        
2011 34865.05 38962.01 33876.42 30611.39 37750.86 30764.14 35867.64 33948.36
          Sep      Oct      Nov      Dec
2010                            56335.97
2011 35359.50 31508.32 31995.30 45813.51

EvoluciĂ³n semanal del consumo. Años 2007-2009

Serie temporal y visualizaciĂ³n

Vamos a representar la serie por semanas. Es decir, para cada submedidor vamos a representar el consumo total semanal. En este caso, la frecuencia serĂ¡ de 7 dĂ­as por año, tenemos una observaciĂ³n para cada dĂ­a de la semana para los años 2007,8 y 2009.

Granularidad_semanas <- Granularidad_horas %>% group_by(year, weekday) %>% 
   summarize(
      Sub_metering_1 = sum (Sub_metering_1),
      Sub_metering_2 = sum (Sub_metering_2),
      Sub_metering_3 = sum (Sub_metering_3),
      Global_active_power = sum (Global_active_power),
      energia2 = sum (energia2)
   )
seriesSemanales <- rbind.data.frame(
   filter(Granularidad_semanas, year== 2007  ),
   filter(Granularidad_semanas, year== 2008  ),
   filter(Granularidad_semanas, year== 2009  )
)

## Creamos un objeto TS
 tsSM1_Semanal<- ts(seriesSemanales$Sub_metering_1,frequency=7, start=c(2007,1))
 tsSM2_Semanal<- ts(seriesSemanales$Sub_metering_2,frequency=7, start=c(2007,1))
 tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3,frequency=7, start=c(2007,1))
## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Semanal, ts.colour = 'red', xlab = "Semana del año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n semanal del consumo de energĂ­a de la cocina \n Años 2007-2009")

autoplot(tsSM2_Semanal, xlab = "Semana del año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n semanal del consumo de energĂ­a de la lavanderĂ­a \n Años 2007-2009")

autoplot(tsSM3_Semanal, ts.colour = 'blue', xlab = "Semana del año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n semanal del consumo de energĂ­a del aire acondicionado y termo elĂ©ctrico \n Años 2007-2009")

Forecasting descomponiendo la serie

DescomposiciĂ³n de la serie y visualizaciĂ³n

Cocina: Submetering 1

Componentes_SM1_SemanalSTL<- tsSM1_Semanal %>%
  stl( t.window=7, s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_SemanalSTL, ts.colour = 'blue')

Se aprecia componente estacional: oscilaciones que se repiten cada año, es decir, cada 7 días. Vemos como la tendencia de consumo es creciente para los años 2007 y 2009 pero sin embargo, en el año 2008 hubo tendencia de consumo decreciente.

plot.ts(tsSM1_Semanal, col = 'gray',main="EvoluciĂ³n del consumo energĂ©tico de cada dĂ­a de la semana de la cocina \n con la tendencia")
lines(Componentes_SM1_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Lavadero: Submetering 2

Componentes_SM2_SemanalSTL<- tsSM2_Semanal %>%stl( t.window=7,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_SemanalSTL, ts.colour = 'blue')

Se observa una tendencia al consumo que es decreciente respecto de los Ăºltimos años, aunque a partir del año 2009 se puede apreciar cierta tendencia creciente del consumo energĂ©tico del lavadero.

Tambien se obserca componente estacional, donde hay un pico de menor consumo energĂ©tico los jueves y domingos, y los sĂ¡bados son lso dĂ­as donde hay un mayor registro del consumo energĂ©tico del lavadero.

plot.ts(tsSM2_Semanal, col = 'gray' , main = 'EvoluciĂ³n del consumo energĂ©tico del lavadero \n por dĂ­a de la semana. Años 2007-2009')
lines(Componentes_SM2_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Aire acondicionado y termo eléctrico: Submetering 3

Componentes_SM3_SemanalSTL<- tsSM3_Semanal %>%stl( t.window=7, 
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_SemanalSTL, ts.colour = 'blue')

Se observa una fuerte componente estacional, con patrones que se repiten para cada año los mismos dĂ­as de la semana. Por ejemplo, los Martes y Viernes son los dĂ­as donde el consumo energĂ©tico es mayor y los sĂ¡bados y MiĂ©rcoles son los dĂ­as donde el consumo energĂ©tico es menor.

La tendencia general del consumo energético es creciente.

plot.ts(tsSM3_Semanal, col = 'gray' , main = 'EvoluciĂ³n del consumo energĂ©tico del AC y termo \n  segĂºn dĂ­a de la semana. Años 2007-2009',)
lines(Componentes_SM3_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Tendencia creciente del consumo energético del aire acondicionado y termo eléctrico

PredicciĂ³n con Holt-Winters (suavizado exponencial) para el año 2010

EvoluciĂ³n semanal del consumo (sin estacionalidad). Años 2007-2009

Cocina
tsSM1_Semanal_detrended <- tsSM1_Semanal - Componentes_SM1_SemanalSTL$time.series[,1]
plot.ts(tsSM1_Semanal_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'EvoluciĂ³n semanal del consumo de la cocina. Años 2007-2009',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM1_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_SemanalSTL$time.series[,3])
qqline(Componentes_SM1_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM1_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM1_SemanalSTL$time.series[, 3]
W = 0.62627, p-value = 3.78e-06

Los resĂ­duos no son solo ruido. AĂºn hay estacionalidad. Volvemos a descomponer la serie desestacionalizada una vez para tratar de eliminar la componente estacional.

## Volvemos a descomponer la serie
plot(stl(tsSM1_Semanal_detrended ,
         s.window = "periodic",robust = TRUE,t.window = 7))

Parece que la componente irregular no sigue ningĂºn patrĂ³n determinado.

Suavizado exponencial de HoltWinters
  • Primera iteracciĂ³n: Solo con \(\alpha\), voy a averiguar en que margen se va a mover la predicciĂ³n
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                                alpha=TRUE
                               # beta=TRUE, 
                               # gamma=TRUE
                                )
summary(tsSM1_HW_Semanal$fitted)
      xhat            level            trend            season        
 Min.   : 47685   Min.   : 67975   Min.   :-900.3   Min.   :-19389.3  
 1st Qu.: 81771   1st Qu.: 80568   1st Qu.:-900.3   1st Qu.: -1906.8  
 Median : 85529   Median : 87483   Median :-900.3   Median :   336.8  
 Mean   : 83958   Mean   : 84859   Mean   :-900.3   Mean   :     0.0  
 3rd Qu.: 92100   3rd Qu.: 89084   3rd Qu.:-900.3   3rd Qu.:  6619.6  
 Max.   :105496   Max.   :101417   Max.   :-900.3   Max.   :  9648.0  
# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronĂ³stico de consumo de la cocina variarĂ¡ con un pronĂ³stico de 47685 vatios hora como mĂ­nimo y 105496 como mĂ¡ximo.

  • Segunda iteracciĂ³n: Añadiendo \(\beta\) podrĂ© averiguar la tendencia del consumo para las predicciones
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                                alpha=TRUE,
                                beta=TRUE
                               # gamma=TRUE
                                )
summary(tsSM1_HW_Semanal$fitted)
      xhat            level            trend              season        
 Min.   : 36880   Min.   : 67975   Min.   :-32119.9   Min.   :-19389.3  
 1st Qu.: 73545   1st Qu.: 80568   1st Qu.: -6256.7   1st Qu.: -1906.8  
 Median : 86848   Median : 87483   Median :  -534.6   Median :   336.8  
 Mean   : 83905   Mean   : 84859   Mean   :  -953.4   Mean   :     0.0  
 3rd Qu.: 94194   3rd Qu.: 89084   3rd Qu.:  2952.5   3rd Qu.:  6619.6  
 Max.   :139838   Max.   :101417   Max.   : 33442.4   Max.   :  9648.0  
# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Parece que los valores extremos tienen gran influencia y se pronosticarĂ¡ un aumento del consumo.

  • Tercera iteracciĂ³n: Si añado \(\gamma\) podrĂ© ver la estacionalidad
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                                alpha=TRUE,
                                beta=TRUE,
                                gamma=TRUE
                                )
summary(tsSM1_HW_Semanal$fitted)
      xhat            level            trend              season        
 Min.   : 36880   Min.   : 67975   Min.   :-32119.9   Min.   :-19389.3  
 1st Qu.: 73545   1st Qu.: 80568   1st Qu.: -6256.7   1st Qu.: -1906.8  
 Median : 86848   Median : 87483   Median :  -534.6   Median :   336.8  
 Mean   : 83905   Mean   : 84859   Mean   :  -953.4   Mean   :     0.0  
 3rd Qu.: 94194   3rd Qu.: 89084   3rd Qu.:  2952.5   3rd Qu.:  6619.6  
 Max.   :139838   Max.   :101417   Max.   : 33442.4   Max.   :  9648.0  
# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 7 observaciones, es decir, la predicciĂ³n del consumo energĂ©tico semanal total para el año 2010.

## HoltWinters forecast & plot
tsSM1_HW_Semanal_for <- forecast(tsSM1_HW_Semanal, h=7)
plot(tsSM1_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha ",
     main="PredicciĂ³n del consumo energĂ©tico semanal de la cocina \n Año 2010" ) 

Vemos como hay un pronĂ³stico de consumo energĂ©tico mucho mayor que para los años anteriores.

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Semanal_forC <- forecast(tsSM1_HW_Semanal, h=7, level = c(95))
## Plot only the forecasted area

plot(tsSM1_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n del consumo energĂ©tico semanal para el año 2010", start(2010))

Lavadero
tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,1]
plot.ts(tsSM2_Semanal_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'EvoluciĂ³n del consumo energĂ©tico segĂºn dĂ­a de la semana \n del lavadero',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM2_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_SemanalSTL$time.series[,3])
qqline(Componentes_SM2_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM2_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM2_SemanalSTL$time.series[, 3]
W = 0.71849, p-value = 4.722e-05

Los resĂ­duos no siguen una distribuciĂ³n normal. La descomposiciĂ³n no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM2_Semanal_detrended ,s.window = "periodic",robust = TRUE,t.window = 7))

Tendencia decreciente del consumo.

Suavizado exponencial de HoltWinters
  • Primera interacciĂ³n: rango de predicciĂ³n
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                               alpha = TRUE
                                )
summary(tsSM2_HW_Semanal$fitted)
      xhat            level            trend           season      
 Min.   :  5231   Min.   : 17137   Min.   :-3959   Min.   :-20239  
 1st Qu.: 65759   1st Qu.: 80171   1st Qu.:-3959   1st Qu.:-13862  
 Median : 95066   Median : 98804   Median :-3959   Median : -7947  
 Mean   : 89219   Mean   : 93179   Mean   :-3959   Mean   :     0  
 3rd Qu.:109099   3rd Qu.:108252   3rd Qu.:-3959   3rd Qu.:  5329  
 Max.   :161540   Max.   :146103   Max.   :-3959   Max.   : 51483  
plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronĂ³stico de consumo energĂ©tico para el lavadero es de 5231 a 161540 vatios hora.

  • Segunda interacciĂ³n: tendencia
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                               alpha = TRUE, beta=TRUE
                                )
summary(tsSM2_HW_Semanal$fitted)
      xhat            level            trend            season      
 Min.   :-67191   Min.   : 17137   Min.   :-76382   Min.   :-20239  
 1st Qu.: 46662   1st Qu.: 80171   1st Qu.:-20120   1st Qu.:-13862  
 Median : 89948   Median : 98804   Median :-10241   Median : -7947  
 Mean   : 91445   Mean   : 93179   Mean   : -1734   Mean   :     0  
 3rd Qu.:142413   3rd Qu.:108252   3rd Qu.: 11518   3rd Qu.:  5329  
 Max.   :262749   Max.   :146103   Max.   :128966   Max.   : 51483  
plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Serie muy afectada por valores extremos, pero se aprecia una tendencia decreciente del consumo energético.

  • Tercera interacciĂ³n
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                               alpha = TRUE,
                               beta=TRUE , gamma=TRUE
                                )
summary(tsSM2_HW_Semanal$fitted)
      xhat            level            trend            season      
 Min.   :-67191   Min.   : 17137   Min.   :-76382   Min.   :-20239  
 1st Qu.: 46662   1st Qu.: 80171   1st Qu.:-20120   1st Qu.:-13862  
 Median : 89948   Median : 98804   Median :-10241   Median : -7947  
 Mean   : 91445   Mean   : 93179   Mean   : -1734   Mean   :     0  
 3rd Qu.:142413   3rd Qu.:108252   3rd Qu.: 11518   3rd Qu.:  5329  
 Max.   :262749   Max.   :146103   Max.   :128966   Max.   : 51483  
plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                             #  alpha = TRUE,
                               beta=TRUE , gamma=TRUE
                                )
summary(tsSM2_HW_Semanal$fitted)
      xhat            level            trend            season        
 Min.   : 39417   Min.   : 60834   Min.   :-11558   Min.   :-23890.7  
 1st Qu.: 58672   1st Qu.: 72658   1st Qu.: -7955   1st Qu.:-14149.9  
 Median : 88846   Median : 85738   Median : -6115   Median : -7378.4  
 Mean   : 85244   Mean   : 90171   Mean   : -3035   Mean   : -1891.6  
 3rd Qu.:101345   3rd Qu.:106734   3rd Qu.:  3060   3rd Qu.:  -584.4  
 Max.   :169618   Max.   :128927   Max.   : 10093   Max.   : 51482.6  
plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Parece que esta configuraciĂ³n se ajuste mejor.

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir el consumo energĂ©tico de los prĂ³ximos 7 dĂ­as, es decir, el consumo energĂ©tico de cada dĂ­a de la semana del año 2010

## HoltWinters forecast & plot
     tsSM2_HW_Semanal_for <- forecast(tsSM2_HW_Semanal, h=7)
plot(tsSM2_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha",
     main="PredicciĂ³n del consumo energĂ©tico del lavadero \ segĂºn dĂ­a de la semana del año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Semanal_forC <- forecast(tsSM2_HW_Semanal, h=7, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha ", main="PredicciĂ³n semanal para el año 2010", start(2010))

AC y termo eléctrico
tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,1]
plot.ts(tsSM3_Semanal_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'EvoluciĂ³n segĂºn dĂ­a de la semana del consumo energĂ©tico del AC y termo \n Serie sin componente estacional',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM3_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_SemanalSTL$time.series[,3])
qqline(Componentes_SM3_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM3_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM3_SemanalSTL$time.series[, 3]
W = 0.5979, p-value = 1.871e-06

Los resĂ­duos no siguen una distribuciĂ³n normal. La descomposiciĂ³n no es buena. Vamos a volver a descomponer la serie.

## Volvemos a descomponer la serie
plot(stl(tsSM3_Semanal_detrended ,s.window = "periodic",robust = TRUE,t.window = 7))

Suavizado exponencial de HoltWinters
  • Primera interacciĂ³n: rango de pronĂ³stico
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                             alpha = TRUE)
summary(tsSM3_HW_Semanal$fitted)
      xhat            level            trend          season        
 Min.   :430345   Min.   :427107   Min.   :3805   Min.   :-14387.8  
 1st Qu.:453104   1st Qu.:446399   1st Qu.:3805   1st Qu.:-10910.1  
 Median :467179   Median :468759   Median :3805   Median :  -566.8  
 Mean   :480676   Mean   :476871   Mean   :3805   Mean   :     0.0  
 3rd Qu.:494414   3rd Qu.:500565   3rd Qu.:3805   3rd Qu.:  8632.7  
 Max.   :566292   Max.   :545292   Max.   :3805   Max.   : 17195.5  
plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PredicciĂ³n del consumo para el termo elĂ©ctrico y aire acondicionado entre 430345 y 566292 vatios hora.

  • Segunda interacciĂ³n: tendencia de consumo
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                             alpha = TRUE,beta=TRUE)
summary(tsSM3_HW_Semanal$fitted)
      xhat            level            trend            season        
 Min.   :419999   Min.   :427107   Min.   :-34437   Min.   :-14387.8  
 1st Qu.:447423   1st Qu.:446399   1st Qu.: -7413   1st Qu.:-10910.1  
 Median :462529   Median :468759   Median :  5721   Median :  -566.8  
 Mean   :485585   Mean   :476871   Mean   :  8714   Mean   :     0.0  
 3rd Qu.:507335   3rd Qu.:500565   3rd Qu.: 27545   3rd Qu.:  8632.7  
 Max.   :600575   Max.   :545292   Max.   : 51106   Max.   : 17195.5  
plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Se observa una tendencia creciente.

  • Tercera iteracciĂ³n
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                             alpha = TRUE,
                             beta=TRUE, gamma = TRUE)
summary(tsSM3_HW_Semanal$fitted)
      xhat            level            trend            season        
 Min.   :419999   Min.   :427107   Min.   :-34437   Min.   :-14387.8  
 1st Qu.:447423   1st Qu.:446399   1st Qu.: -7413   1st Qu.:-10910.1  
 Median :462529   Median :468759   Median :  5721   Median :  -566.8  
 Mean   :485585   Mean   :476871   Mean   :  8714   Mean   :     0.0  
 3rd Qu.:507335   3rd Qu.:500565   3rd Qu.: 27545   3rd Qu.:  8632.7  
 Max.   :600575   Max.   :545292   Max.   : 51106   Max.   : 17195.5  
plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 7 observaciones, es decir, el consumo energĂ©tico de cada dĂ­a de la semana para el año 2010.

## HoltWinters forecast & plot
tsSM3_HW_Semanal_for <- forecast(tsSM3_HW_Semanal, h=7)
plot(tsSM3_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha",
     main="PredicciĂ³n del consumo energĂ©tico semanal \n del AC y termo elĂ©ctrico. Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Semanal_forC <- forecast(tsSM3_HW_Semanal, h=7, level=c(95))
## Plot only the forecasted area
plot(tsSM3_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n semanal del consumo energĂ©tico para el año 2010", start(2010), ylim=c(0,450000))

ComparaciĂ³n de los coeficientes y resultados de cada predicciĂ³n. Datos sin estacionalidad

En las grĂ¡ficas de predicciĂ³n, hemos observado que el ajuste no mantiene la tendencia de consumo energĂ©tico de la serie.

DatosRealesSeriesSemanales<-Granularidad_semanas  %>% filter(year==2010)
RMSE_SM1_GranSemanalFor<-accuracy(tsSM1_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
RMSE_SM2_GranSemanalFor<-accuracy(tsSM2_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
RMSE_SM3_GranSemanalFor<-accuracy(tsSM3_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_3)[4]


ResEvolSemanal2_SinEstacionalidadConTendencia <-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranSemanalFor,
               RMSE_SM2_GranSemanalFor,
               RMSE_SM3_GranSemanalFor)
)
ResEvolSemanal2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 28600.80
Lavadero 51998.21
AC y TermoE 218301.40

Los errores son demasiado altos, las predicciones no serĂ¡n del todo fiables y por tanto no podemos ofrecerlas al cliente.

AnĂ¡lisis de series temporales con metodologĂ­a Box-Jenkis

Paso 1. Lectura y representaciĂ³n grĂ¡fica de los datos

# tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3,         frequency=7, start=c(2007,1))
autoplot(tsSM3_Semanal, ts.colour = 'blue', 
         xlab = "Semana del año",
         ylab = "EnergĂ­a (Varios-hora)", 
         main = "EvoluciĂ³n semanal del consumo de energĂ­a del aire acondicionado y termo elĂ©ctrico \n Años 2007-2009")

Paso 2. Transformaciones para que la varianza sea estable en el tiempo

library(TSA)
bc=BoxCox.ar(y=tsSM3_Semanal);

bc$mle
[1] 2
bc$ci # El 1 no estĂ¡ en el intervalo, hay que transformar los datos
[1] -1.1  2.0
# El 0 no estĂ¡ en el intervalo, no es transformaciĂ³n logarĂ­tmica

Usando la familia de transformaciones BoxCox, nos sugiere un valor de lambda = 2. El valor 1 (no transformar) estĂ¡ en el intervalo, por lo que no hay que hacer una transformaciĂ³n de los datos, podemos asumir que la varianza es constante en el tiempo

# TransftsSM3_Semanal<- tsSM3_Semanal^2

Paso 3. Transformaciones para que la media sea estable en el tiempo

acf(tsSM3_Semanal, main="FAS del consumo de EnergĂ­a del AC y Termo", lag.max=50)

ndiffs( tsSM3_Semanal) # Necesito una diferencia regular
[1] 1
nsdiffs(tsSM3_Semanal)
[1] 0

Diferencia estacional:

tsSM3_SemanalDif=diff(tsSM3_Semanal, lag=1, diff=1)
acf(tsSM3_SemanalDif, main="FAS de primera diferencia regular del consumo energético", lag=50)

ndiffs(tsSM3_SemanalDif) 
[1] 0

Parece la FAS de un Ruido Blanco.

Paso 4. Contrastar la estacionariedad

Test de raiz unitaria.

library(tseries)
adf.test(tsSM3_SemanalDif)

    Augmented Dickey-Fuller Test

data:  tsSM3_SemanalDif
Dickey-Fuller = -6.2783, Lag order = 2, p-value = 0.01
alternative hypothesis: stationary

No existen evidencias significativas para afirmar que los datos no sean estacionarios. La media es estable en el tiempo.

autoplot(tsSM3_SemanalDif, ts.colour = 'blue', 
         xlab = "Semana del año",
         ylab = "EnergĂ­a (Varios-hora)", 
         main = "EvoluciĂ³n semanal del consumo de energĂ­a del aire acondicionado y termo elĂ©ctrico \n Años 2007-2009")

Paso 5. Identificar la estructura ARIMA

acf(tsSM3_SemanalDif, main="FAS tras una diferencia regular", lag=50)

Ruido blanco

pacf(tsSM3_SemanalDif, main="FAP tras una diferencia regular", lag=50)

AR(1).

Pasos 6 y 7. EstimaciĂ³n de parĂ¡metros y diagnĂ³stico

modelo ARIMA(1,1,0)xARIMA(0,0,0)

ajuste1=arima(tsSM3_Semanal, order=c(1,1,0), seasonal=list(order=c(0,0,0), period=7)); ajuste1

Call:
arima(x = tsSM3_Semanal, order = c(1, 1, 0), seasonal = list(order = c(0, 0, 
    0), period = 7))

Coefficients:
          ar1
      -0.4254
s.e.   0.2312

sigma^2 estimated as 2.967e+09:  log likelihood = -246.59,  aic = 495.17
checkresiduals(ajuste1)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)
Q* = 14.368, df = 3, p-value = 0.002445

Model df: 1.   Total lags used: 4

No pasa la diagnosis: p-value = 0.002445. El modelo no es vĂ¡lido.

ajuste2=auto.arima(tsSM3_Semanal, d = 1,D = 0); ajuste2
Series: tsSM3_Semanal 
ARIMA(3,1,0) with drift 

Coefficients:
          ar1      ar2      ar3     drift
      -1.2341  -1.1548  -0.6579  6092.903
s.e.   0.2171   0.1892   0.2114  1912.826

sigma^2 = 1.26e+09:  log likelihood = -237.16
AIC=484.31   AICc=488.6   BIC=489.29

Modelo propuesto: ARIMA(3,1,0)

aicAutoArima = 496.21 AICModelo1=484.31

checkresiduals(ajuste2)


    Ljung-Box test

data:  Residuals from ARIMA(3,1,0) with drift
Q* = 2.8112, df = 3, p-value = 0.4217

Model df: 4.   Total lags used: 7

El modelo pasa la diagnosis, p-valor=0.422, por tanto es vĂ¡lido. No existen evidencias significativas en contra de que los resĂ­duos sigan un proceso de ruido blanco.

Nos quedamos con el modelo 2.

Paso 8. PredicciĂ³n de los resultados

pred=forecast(ajuste2,h = 7)# predicciĂ³n del año 2010
plot(tsSM3_Semanal, xlim=c(2007, 2011), ylim=c(350000,600000))
lines(pred$mean, col="red")

# plot((forecast(ajuste.bueno,h=7))) #muy similar al anterior

Vamos a calcular el error: 4.9418866^{5}, que es mucho mayor que por el otro método.

PronĂ³stico de consumo para los prĂ³ximos 12 meses

Lo haremos con el modelo conseguido primero, no con la metodologĂ­a de Box-Jenkis, que es peor.

Cocina

# DefiniciĂ³n de serie temporal

Granularidad_semanas <- Granularidad_horas %>% group_by(year, weekday) %>% 
   summarize(
      Sub_metering_1 = sum (Sub_metering_1),
      Sub_metering_2 = sum (Sub_metering_2),
      Sub_metering_3 = sum (Sub_metering_3),
      Global_active_power = sum (Global_active_power),
      energia2 = sum (energia2)
   )
seriesSemanalesCompletas <- Granularidad_semanas

## Creamos un objeto TS
 tsSM1_Semanal<- ts(seriesSemanalesCompletas$Sub_metering_1,frequency=7, start=c(2007,1))
 tsSM2_Semanal<- ts(seriesSemanalesCompletas$Sub_metering_2,frequency=7, start=c(2007,1))
 tsSM3_Semanal<- ts(seriesSemanalesCompletas$Sub_metering_3,frequency=7, start=c(2007,1))
# DescomposiciĂ³n de la serie
Componentes_SM1_SemanalSTL<-  tsSM1_Semanal %>%
  stl( t.window=7,s.window="periodic",  robust=TRUE)
# EliminaciĂ³n de la estacionalidad
tsSM1_Semanal_detrended <- tsSM1_Semanal- Componentes_SM1_SemanalSTL$time.series[,1]
# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended , 
                                seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)
# PronĂ³stico para los siguientes 13 meses 
PrediccionFinalCocinaSemanal <- forecast(tsSM1_HW_Semanal, h=7)
## Plot only the forecasted area
plot(PrediccionFinalCocinaSemanal, 
     ylab= "Vatios-Hora", xlab="Fecha", 
     main="PredicciĂ³n el consumo energĂ©tico de la cocina \n segĂºn dĂ­a de para año 2011", start(2010,12), ylim = c(100000,200000))

Las predicciones son las siguientes:

PrediccionFinalCocinaSemanal[["mean"]] 
Time Series:
Start = c(2012, 1) 
End = c(2012, 7) 
Frequency = 7 
[1] 173867.6 190256.8 182292.0 184206.4 199466.0 133717.8 125668.2

SegĂºn el pronĂ³stico, el dĂ­a de la semana que mĂ¡s consumo energĂ©tico total de energĂ­a habrĂ¡, serĂ¡ el Viernes. Y por el contrario, el dĂ­a de la semana en que se registrarĂ¡ el menor consumo energĂ©tico es el Domingo.

Lavadero

# DescomposiciĂ³n de la serie
Componentes_SM2_SemanalSTL<-  
  tsSM2_Semanal %>%
  stl( t.window=7,s.window="periodic",  robust=TRUE)
# EliminaciĂ³n de la estacionalidad
tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,1]
# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended , 
                                seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE
 #  alpha=TRUE
)
# PronĂ³stico para los siguientes 13 meses 
PrediccionFinalLavaderoSemanal <- forecast(tsSM2_HW_Semanal, h=7)
## Plot only the forecasted area
plot(PrediccionFinalLavaderoSemanal, ylab= "Vatios-Hora", xlab="Fecha", main="PredicciĂ³n el consumo energĂ©tico del lavadero \n  SegĂºn el dĂ­a de la semana para los prĂ³ximos 12 meses", start(2010,12), ylim = c(50000,450000))

Las predicciones son las siguientes:

PrediccionFinalLavaderoSemanal[["mean"]] 
Time Series:
Start = c(2012, 1) 
End = c(2012, 7) 
Frequency = 7 
[1]  67555.2 134504.6 295918.7 390770.7 420871.3 451580.6 452544.9

Aire acondicionado y termo eléctrico

# DescomposiciĂ³n de la serie
Componentes_SM3_SemanalSTL<-  tsSM3_Semanal %>%
  stl( t.window=7,s.window="periodic",  robust=TRUE)
# EliminaciĂ³n de la estacionalidad
tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,1]
# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended , 
                                seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)
# PronĂ³stico para los siguientes 13 meses 
PrediccionFinalAcTSemanal <- forecast(tsSM3_HW_Semanal, h=13)
## Plot only the forecasted area
plot(PrediccionFinalAcTSemanal, ylab= "Vatios-Hora", xlab="Fecha", main="PredicciĂ³n el consumo energĂ©tico del aire y termo \n para cada dĂ­a de la semana", start(2010,12), ylim = c(350000,850000))

Las predicciones son las siguientes:

PrediccionFinalLavaderoSemanal[["mean"]] 
Time Series:
Start = c(2012, 1) 
End = c(2012, 7) 
Frequency = 7 
[1]  67555.2 134504.6 295918.7 390770.7 420871.3 451580.6 452544.9

EvoluciĂ³n del consumo energĂ©tico por semana del año.

Serie temporal y visualizaciĂ³n

Vamos a representar la serie por meses para cada submedidor y tambien para la energĂ­a global co granularidad mensual. En este caso, la frecuencia serĂ¡ 12, tenemos una observaciĂ³n por mes para los años 2007,8 y 2009.

Granularidad_semanas <- Granularidad_horas %>% group_by(year, week) %>% 
   summarize(
      Sub_metering_1 = sum (Sub_metering_1),
      Sub_metering_2 = sum (Sub_metering_2),
      Sub_metering_3 = sum (Sub_metering_3),
      Global_active_power = sum (Global_active_power),
      energia2 = sum (energia2)
   )
seriesSemanales <- rbind.data.frame(
   filter(Granularidad_semanas, year== 2007  ),
   filter(Granularidad_semanas, year== 2008  ),
   filter(Granularidad_semanas, year== 2009  )
)

## Creamos un objeto TS
   tsSM1_Semanal<- ts(seriesSemanales$Sub_metering_1,         frequency=53, start=c(2007,1))
   tsSM2_Semanal<- ts(seriesSemanales$Sub_metering_2,         frequency=53, start=c(2007,1))
   tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3,         frequency=53, start=c(2007,1))
tsSM_GAP_Semanal<- ts(seriesSemanales$Global_active_power,    frequency=53, start=c(2007,1))
## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Semanal, ts.colour = 'red', xlab = "Semana del año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n semanal del consumo de energĂ­a de la cocina \n Años 2007-2009")

autoplot(tsSM2_Semanal, xlab = "Semana del año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n semanal del consumo de energĂ­a de la lavanderĂ­a \n Años 2007-2009")

autoplot(tsSM3_Semanal, ts.colour = 'blue', xlab = "Semana del año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n semanal del consumo de energĂ­a del aire acondicionado y termo elĂ©ctrico \n Años 2007-2009")

autoplot(tsSM_GAP_Semanal, ts.colour = 'green', xlab = "Semana del año", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n semanal del consumo de energĂ­a \n Años 2007-2009")

AC y Termo: podrĂ­a apreciarse cierta estacionalidad EnergĂ­a total: muelle

PredicciĂ³n del año 2010 antes de descomponer la serie

Vamos a plicar regresiĂ³n lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrĂ¡tico medio.

SM1. Cocina

Modelo

fitSM1_Semanal <- tslm(tsSM1_Semanal ~ trend + season) 
summary(fitSM1_Semanal)

Call:
tslm(formula = tsSM1_Semanal ~ trend + season)

Residuals:
     Min       1Q   Median       3Q      Max 
-10942.7  -1970.9   -182.3   2205.3  12203.0 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 10835.442   2693.122   4.023 0.000108 ***
trend          -8.872      8.426  -1.053 0.294745    
season2      6872.206   3753.914   1.831 0.069985 .  
season3      6643.745   3753.942   1.770 0.079662 .  
season4      4273.617   3753.989   1.138 0.257536    
season5      3290.489   3754.056   0.877 0.382751    
season6      4326.029   3754.141   1.152 0.251800    
season7       791.901   3754.245   0.211 0.833347    
season8      -313.227   3754.368  -0.083 0.933669    
season9     -3984.688   3754.509  -1.061 0.290985    
season10     6229.851   3754.670   1.659 0.100054    
season11     2045.390   3754.850   0.545 0.587092    
season12     5559.929   3755.048   1.481 0.141693    
season13     3915.802   3755.266   1.043 0.299459    
season14     1586.341   3755.502   0.422 0.673594    
season15     3100.547   3755.757   0.826 0.410934    
 [ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4598 on 105 degrees of freedom
Multiple R-squared:  0.4836,    Adjusted R-squared:  0.223 
F-statistic: 1.855 on 53 and 105 DF,  p-value: 0.003663

Resultados:

  • Valor del coeficiente de determinaciĂ³n: 0.4836. Es un valor bajo, el modelo Ăºnicamente explica el 48% de la variabilidad total.

  • Hay un coeficiente significativamente no nulo.

PredicciĂ³n IC

PredicciĂ³n para los prĂ³ximos 20 dĂ­as

forecastfitSM1_Semanal <- forecast(fitSM1_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="PredicciĂ³n del consumo energĂ©tico semanal total del año 2010 \n  a travĂ©s de un modelo de regresiĂ³n")

SM2. Lavadero

Modelo

fitSM2_Semanal <- tslm(tsSM2_Semanal ~ trend + season) 
summary(fitSM2_Semanal)

Call:
tslm(formula = tsSM2_Semanal ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-12710  -2492    262   2070   9448 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  19974.287   2851.564   7.005 2.45e-10 ***
trend          -46.079      8.921  -5.165 1.15e-06 ***
season2      -4864.254   3974.765  -1.224 0.223773    
season3        714.159   3974.795   0.180 0.857756    
season4       -986.429   3974.845  -0.248 0.804489    
season5      -1733.682   3974.915  -0.436 0.663619    
season6       -145.270   3975.005  -0.037 0.970917    
season7      -4670.524   3975.115  -1.175 0.242676    
season8      -1727.111   3975.245  -0.434 0.664841    
season9      -5223.698   3975.395  -1.314 0.191707    
season10      3920.048   3975.565   0.986 0.326383    
season11      1273.127   3975.756   0.320 0.749436    
season12     -1349.793   3975.966  -0.339 0.734920    
season13      -436.381   3976.196  -0.110 0.912818    
season14     -2924.635   3976.446  -0.735 0.463681    
season15     -2526.222   3976.716  -0.635 0.526645    
 [ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4868 on 105 degrees of freedom
Multiple R-squared:  0.5326,    Adjusted R-squared:  0.2966 
F-statistic: 2.257 on 53 and 105 DF,  p-value: 0.0001984

Resultados:

  • Valor del coeficiente de determinaciĂ³n: 0.5326. Es un valor aceptable, el modelo Ăºnicamente explica el 53% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

PredicciĂ³n IC

PredicciĂ³n para los prĂ³ximos 12 meses (año 2010)

forecastfitSM2_Semanal <- forecast(fitSM2_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="PredicciĂ³n del consumo energĂ©tico del año 2010 \n  a travĂ©s de un modelo de regresiĂ³n")

Se prevee una bajada del consumo

SM3. AC y termo

Modelo

fitSM3_Semanal <- tslm(tsSM3_Semanal ~ trend + season) 
summary(fitSM3_Semanal)

Call:
tslm(formula = tsSM3_Semanal ~ trend + season)

Residuals:
   Min     1Q Median     3Q    Max 
-35933  -5644    -17   6305  26151 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  56824.99    7005.03   8.112 9.94e-13 ***
trend           95.69      21.92   4.366 2.97e-05 ***
season2      13899.65    9764.23   1.424  0.15755    
season3      18434.30    9764.31   1.888  0.06180 .  
season4       9370.94    9764.43   0.960  0.33941    
season5      20968.59    9764.60   2.147  0.03406 *  
season6      14556.57    9764.82   1.491  0.13903    
season7       8915.22    9765.09   0.913  0.36335    
season8      -2728.46    9765.41  -0.279  0.78049    
season9     -25352.82    9765.78  -2.596  0.01078 *  
season10      2945.17    9766.20   0.302  0.76358    
season11      8872.81    9766.67   0.908  0.36571    
season12      8371.46    9767.18   0.857  0.39334    
season13     13322.78    9767.75   1.364  0.17550    
season14      9586.09    9768.36   0.981  0.32868    
season15     -1125.59    9769.03  -0.115  0.90849    
 [ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 11960 on 105 degrees of freedom
Multiple R-squared:  0.7153,    Adjusted R-squared:  0.5715 
F-statistic: 4.977 on 53 and 105 DF,  p-value: 1.341e-12

Resultados:

  • Valor del coeficiente de determinaciĂ³n: 0.7153. Es un valor que podrĂ­a ser aceptable, el modelo Ăºnicamente explica el 71% de la variabilidad total.

  • Hay tres coeficientes significativamente no nulos.

PredicciĂ³n IC

PredicciĂ³n para los prĂ³ximos 12 meses (año 2010)

forecastfitSM3_Semanal <- forecast(fitSM3_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Semanal, ylab= "Watt-Hours", xlab="Fecha",
     main="PredicciĂ³n del consumo energĂ©tico del año 2010 \n  a travĂ©s de un modelo de regresiĂ³n")

Se prevee una bajada del consumo

ComparaciĂ³n de los coeficientes y resultados de cada modelo

###  DatosRealesSeriesSemanales <- Granularidad_semanas %>% filter(year == 2010 ###  )
###  
###  RMSE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
###    R2_SM1_GranSemanal<-0.975
###  MASE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_1)[12]
###  RMSE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
###    R2_SM2_GranSemanal<-0.7717
###  MASE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_2)[12]
###  RMSE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_3)[4]
###    R2_SM3_GranSemanal<-0.7601
###  MASE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal  ###  ,DatosRealesSeriesSemanales$Sub_metering_3)[12]
###  
###  ResEvolSemanal<-
###  cbind.data.frame(
###        Submedidor=c("Cocina","Lavadero","AC y TermoE"),
###        RMSE = c(RMSE_SM1_GranSemanal,
###                 RMSE_SM2_GranSemanal,
###                 RMSE_SM3_GranSemanal),
###       # MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMen###  sual),
###        R2 = c(R2_SM1_GranSemanal,
###               R2_SM2_GranSemanal,
###               R2_SM3_GranSemanal))
###  ResEvolSemanal %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = ###  "striped")
###  
###  # MASE: Mean absolute error
###  # RMSE: root mean scuare error. Diferencia entre los valores predichos y ###  los observados.

El RMSE del modelo con mejor coeficiente de determinaciĂ³n es el mĂ¡s alto. Pero no hay que olvidar que el consumo elĂ©ctrico para el submedidor correspondiente al AC y termo elĂ©ctrico es el que mĂ¡s energĂ­a consume de todos con una diferencia muy grande.

Forecasting descomponiendo la serie (asĂ­ ok)

Descomponer: quitar tendencia

DescomposiciĂ³n de la serie y visualizaciĂ³n

Submetering 1

 tsSM1_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal              trend             remainder         
 Min.   :-11380.560   Min.   : 9373.389   Min.   :-17306.164  
 1st Qu.: -1525.527   1st Qu.:10411.060   1st Qu.: -1565.107  
 Median :   605.971   Median :11525.838   Median :    50.446  
 Mean   :     0.000   Mean   :11285.767   Mean   :   161.290  
 3rd Qu.:  2468.518   3rd Qu.:12047.600   3rd Qu.:  1583.123  
 Max.   : 12092.727   Max.   :12804.151   Max.   : 18204.879  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     3994         1637      3148          6386 
   %  62.5         25.6      49.3         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8303  0.9452  0.8244  0.9824  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 1591 81 53
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 160 9 6
 $ inner: int 1
 $ outer: int 15
#tsSM1_Semanal %>%
#  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
#  autoplot()

Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

Componentes_SM1_SemanalSTL<- tsSM1_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM1_Semanal, col = 'gray')
lines(Componentes_SM1_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Submetering 2

 tsSM2_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE) %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal              trend            remainder         
 Min.   :-11082.099   Min.   :10972.07   Min.   :-17744.999  
 1st Qu.: -2312.887   1st Qu.:11343.85   1st Qu.: -1509.064  
 Median :   306.932   Median :12813.79   Median :   128.053  
 Mean   :     0.000   Mean   :13181.15   Mean   :   101.019  
 3rd Qu.:  2737.120   3rd Qu.:15219.46   3rd Qu.:  1570.218  
 Max.   :  9565.068   Max.   :15803.51   Max.   : 14301.230  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     5050         3876      3079          7916 
   %  63.8         49.0      38.9         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.7600  0.9452  0.7848  0.9927  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 1591 81 53
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 160 9 6
 $ inner: int 1
 $ outer: int 15
#tsSM2_Semanal %>%
#  stl(t.window=13, s.window="periodic", robust=TRUE) %>%
#  autoplot()
Componentes_SM2_SemanalSTL<- tsSM2_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM2_Semanal, col = 'gray')
lines(Componentes_SM2_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

Submetering 3

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente

tsSM3_Semanal %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend            remainder        
 Min.   :-47211.03   Min.   :57753.65   Min.   :-56145.14  
 1st Qu.:-10494.71   1st Qu.:58324.34   1st Qu.: -4929.51  
 Median :  4528.45   Median :63141.11   Median :  -321.52  
 Mean   :     0.00   Mean   :62123.83   Mean   :  -707.60  
 3rd Qu.: 11265.72   3rd Qu.:63648.44   3rd Qu.:  3997.01  
 Max.   : 26621.30   Max.   :68148.68   Max.   : 39973.53  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     21760         5324      8927         23820
   %  91.4         22.4      37.5         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8598  0.9452  0.8486  0.9898  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 1591 81 53
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 160 9 6
 $ inner: int 1
 $ outer: int 15
Componentes_SM3_SemanalSTL<- tsSM3_Semanal %>%stl( #t.window=13,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_SemanalSTL, ts.colour = 'blue')

plot.ts(tsSM3_Semanal, col = 'gray')
lines(Componentes_SM3_SemanalSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Tendencia creciente del consumo energético

PredicciĂ³n con Holt-Winters (suavizado exponencial) para el año 2010

Sin estacionalidad

Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.

Tendencia de predicciĂ³n (crecimiento o decrecimiento del consumo) y predicciĂ³n teniendo en cuenta las componentes

EvoluciĂ³n mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina
#tsSM1_Mensual_Adjusted <- tsSM1_Mensual - Componentes_SM1_Mensual$seasonal
#autoplot(tsSM1_Mensual_Adjusted)
tsSM1_Semanal_detrended <- tsSM1_Semanal - Componentes_SM1_SemanalSTL$time.series[,2]
plot.ts(tsSM1_Semanal_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM1_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_SemanalSTL$time.series[,3])
qqline(Componentes_SM1_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM1_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM1_SemanalSTL$time.series[, 3]
W = 0.83235, p-value = 3.231e-12

Los resĂ­duos no son solo ruido. AĂºn hay estacionalidad

Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
                              #  tsSM1_Mensual_Adjusted, 
  beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 53 observaciones, es decir, el consumo energĂ©tico semanal de cada submedidor para las 53 semanas del año 2010.

## HoltWinters forecast & plot
tsSM1_HW_Semanal_for <- forecast(tsSM1_HW_Semanal, h=53)
plot(tsSM1_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="PredicciĂ³n del consumo energĂ©tico semanal de la cocina \n Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Semanal_forC <- forecast(tsSM1_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n del consumo energĂ©tico semanal para el año 2010", start(2010))

Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM2_Mensual_Adjusted <- tsSM2_Mensual - Componentes_SM2_Mensual$seasonal
# autoplot(tsSM2_Mensual_Adjusted)
tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,2]
plot.ts(tsSM2_Semanal_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM2_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_SemanalSTL$time.series[,3])
qqline(Componentes_SM2_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM2_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM2_SemanalSTL$time.series[, 3]
W = 0.88071, p-value = 5.464e-10

Los resĂ­duos no siguen una distribuciĂ³n normal. La descomposiciĂ³n no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM2_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
                                # tsSM2_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 7 observaciones.

## HoltWinters forecast & plot
     tsSM2_HW_Semanal_for <- forecast(tsSM2_HW_Semanal, h=53)
plot(tsSM2_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n semanal del consumo energĂ©tico del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Semanal_forC <- forecast(tsSM2_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n semanal para el año 2010", start(2010))

AC y termo eléctrico
# ## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM3_Mensual_Adjusted <- tsSM3_Mensual - Componentes_SM3_Mensual$seasonal
# autoplot(tsSM3_Mensual_Adjusted)
tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,2]
plot.ts(tsSM3_Semanal_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'Seasonal variability of energy comsumption',
        cex.main = 0.85)

par(mfrow = c(1,2))
plot(Componentes_SM3_SemanalSTL$time.series[,3], 
     col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_SemanalSTL$time.series[,3])
qqline(Componentes_SM3_SemanalSTL$time.series[,3])

shapiro.test(Componentes_SM3_SemanalSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM3_SemanalSTL$time.series[, 3]
W = 0.83537, p-value = 4.311e-12

Los resĂ­duos no siguen una distribuciĂ³n normal. La descomposiciĂ³n no es buena.

## Volvemos a descomponer la serie
plot(stl(tsSM3_Semanal_detrended ,s.window = "periodic",robust = TRUE))

La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
                              # tsSM3_Mensual_Adjusted,
                                beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Semanal,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Semanal_for <- forecast(tsSM3_HW_Semanal, h=53)
plot(tsSM3_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n del consumo energĂ©tico semanal del AC y termo elĂ©ctrico \n Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Semanal_forC <- forecast(tsSM3_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n semanal del consumo energĂ©tico para el año 2010", start(2010))

ComparaciĂ³n de los coeficientes y resultados de cada predicciĂ³n. Datos sin estacionalidad

En las grĂ¡ficas de predicciĂ³n, hemos observado que el ajuste no mantiene la tendencia de consumo energĂ©tico de la serie.

RMSE_SM1_GranSemanalFor<-accuracy(tsSM1_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
RMSE_SM2_GranSemanalFor<-accuracy(tsSM2_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
RMSE_SM3_GranSemanalFor<-accuracy(tsSM3_HW_Semanal_forC   ,DatosRealesSeriesSemanales$Sub_metering_3)[4]


ResEvolSemanal2_SinEstacionalidadConTendencia <-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranSemanalFor,
               RMSE_SM2_GranSemanalFor,
               RMSE_SM3_GranSemanalFor)
)
ResEvolSemanal2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 64689.54
Lavadero 91145.26
AC y TermoE 545792.05

EvoluciĂ³n mensual del consumo (con estacionalidad). Años 2007-2009

Vamos a ver una predicciĂ³n del consumo enegĂ©tico de cada submedidor con el mĂ©todo de HW. En este caso, no eliminaremos la estacionalidad, para ver asĂ­ si se mantiene la tendencia de consumo.

Cocina
tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)

Parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM1_Mensual_AdjustedII))

La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
     main="PredicciĂ³n del consumo energĂ©tico de la cocina \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n para el año 2010", start(2010))

Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)

No arece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM2_Mensual_AdjustedII)) 

TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=53)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n del consumo energĂ©tico del lavadero \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n para el año 2010", start(2010))

AC y termo eléctrico
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)

No parece un muelle.

## Volvemos a descomponer la serie
plot(decompose(tsSM3_Mensual_AdjustedII))

Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_MensualII,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=53)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n del consumo energĂ©tico del AC y termo elĂ©ctrico \ Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n para el año 2010", start(2010))

ComparaciĂ³n de los coeficientes y resultados de cada predicciĂ³n. Datos con estacionalidad
RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII   ,DatosRealesSeriesMensuales$Sub_metering_3)[4]


ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_GranMensualForII,
               RMSE_SM2_GranMensualForII,
               RMSE_SM3_GranMensualForII)
   
    
)

ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 18462.43
Lavadero 36358.95
AC y TermoE 433409.63

ComparaciĂ³n de resultados. PredicciĂ³n antes y despues de eliminar estacionalidad

Comparacion<-cbind.data.frame(
  ResEvolMensual2_SinEstacionalidadConTendencia,
  ResEvolMensual2SinTendenciaConEstacionalidad[,2])

colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE Con Estacionalidad RMSE Sin estacionalidad
Cocina 14278.84 18462.43
Lavadero 44294.27 36358.95
AC y TermoE 325349.90 433409.63

Resultados mĂ¡s fiables y con menor error trĂ¡s haber eliminado la estacionalidad de la serie. TambiĂ©n asĂ­ podemos ver la tendencia.

#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
#             grobs = list(p1,p2),
#           nrow = 2,
#             top="Sin estacionalidad",bottom="Con estacionalidad"
#             )

EvoluciĂ³n diaria del consumo por meses. Mes de Enero

Serie temporal y visualizaciĂ³n

En este caso, crearemos una serie para cada mes, ya que la frecuencia de los datos no serĂ¡ la misma, al tener cada mes un nĂºmero de dĂ­as diferente.

par(mfrow=c(2,2))
seriesDiariasEnero <- Granularidad_dias %>% filter(year != 2006 & year != 2010 & month == 1)
## Creamos un objeto TS
   tsSM1_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_1, frequency=31, start=c(2007,1))
   
   tsSM2_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_2, frequency=31, start=c(2007,1))
   
   tsSM3_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_3, frequency=31, start=c(2007,1))
## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_DiariaEnero, ts.colour = 'red', xlab = "DĂ­a", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n diaria del consumo de energĂ­a de la cocina en el mes de Enero \n Años 2007-2009")

autoplot(tsSM2_DiariaEnero, xlab = "DĂ­a", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n diaria del consumo de energĂ­a de la lavanderĂ­a en el mes de Enero \n Años 2007-2009")

autoplot(tsSM3_DiariaEnero, ts.colour = 'blue', xlab = "DĂ­a", ylab = "EnergĂ­a (Varios-hora)", main = "EvoluciĂ³n diaria del consumo energĂ©tico del aire acondicionado y termo en el mes de Enero \n Años 2007-2009")

Parece que no existen muchos patrones que nos podrĂ­an ayudar a hacer predicciones, no serĂ­an muy fiables. (grĂ¡fica muelle)

Forecasting descomponiendo la serie (asĂ­ ok)

DescomposiciĂ³n de la serie y visualizaciĂ³n

Cocina: Submetering 1

Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.

 tsSM1_DiariaEnero %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend             remainder        
 Min.   :-1324.873   Min.   : 871.8631   Min.   :-6121.092  
 1st Qu.: -761.466   1st Qu.:1309.3572   1st Qu.: -491.953  
 Median : -158.424   Median :1491.6212   Median :  220.447  
 Mean   :    0.000   Mean   :1522.7948   Mean   :  550.818  
 3rd Qu.:  337.611   3rd Qu.:1775.0069   3rd Qu.:  732.488  
 Max.   : 5397.352   Max.   :2077.1908   Max.   : 9934.040  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
     1099.1        465.6    1224.4        1692.0
   %  65.0         27.5      72.4         100.0 

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8124  0.9452  0.7973  0.9880  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 931 47 31
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 94 5 4
 $ inner: int 1
 $ outer: int 15

La componente irregular tiene bastante peso, al igual que la estacional.

Componentes_SM1_DiariaEneroSTL<-
  tsSM1_DiariaEnero %>%stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_DiariaEneroSTL, ts.colour = 'blue')

La tendencia del consumo energético de la cocina para los meses de Enero es creciente. Se observa una fuerte componenete estacional.

plot.ts( tsSM1_DiariaEnero, col = 'gray',main="EvoluciĂ³n del consumo energĂ©tico del lavadero \n  de los meses de Enero")
lines(Componentes_SM1_DiariaEneroSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Lavadero: Submetering 2

 tsSM2_DiariaEnero %>%
  stl( #t.window=13,
      s.window="periodic", 
      robust=TRUE) %>% summary()
 Call:
 stl(x = ., s.window = "periodic", robust = TRUE)

 Time.series components:
    seasonal             trend            remainder        
 Min.   :-1788.921   Min.   :1491.391   Min.   :-5599.893  
 1st Qu.:-1294.941   1st Qu.:1617.887   1st Qu.: -685.829  
 Median : -782.743   Median :1866.286   Median : -122.576  
 Mean   :    0.000   Mean   :1939.768   Mean   :  373.264  
 3rd Qu.:  718.901   3rd Qu.:2339.704   3rd Qu.: 1313.902  
 Max.   : 4203.925   Max.   :2437.168   Max.   : 6567.217  
 IQR:
     STL.seasonal STL.trend STL.remainder data  
     2013.8        721.8    1999.7        3494.0
   %  57.6         20.7      57.2         100.0 

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8067  0.9452  0.7789  0.9886  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 931 47 31
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 94 5 4
 $ inner: int 1
 $ outer: int 15
Componentes_SM2_DiariaEneroSTL<- tsSM2_DiariaEnero %>%
  stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_DiariaEneroSTL, ts.colour = 'blue')

No hay una Ăºnica tendecia, los primeros dĂ­as del mes de Enero del año 2007 se aprecia un aumento del consumo energĂ©tico, al ogual que ocurre desde la segunda mitad del mes de enero del año 2008 hasta finales de 2009. Sin embargo, hay un descenso del consumo desde la segunda mitad del año 2007 hasta la misma fecha del año 2008

plot.ts(tsSM2_DiariaEnero, col = 'gray', main="EvoluciĂ³n del consumo energĂ©tico \n  del lavadero el mes de Enero")
lines(Componentes_SM2_DiariaEneroSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

Aire acondicionado y termo eléctrico: Submetering 3

tsSM3_DiariaEnero %>%
  stl( t.window=31,
      s.window="periodic", 
      robust=TRUE)  %>% summary()
 Call:
 stl(x = ., s.window = "periodic", t.window = 31, robust = TRUE)

 Time.series components:
    seasonal             trend             remainder        
 Min.   :-3370.606   Min.   : 7519.501   Min.   :-9218.843  
 1st Qu.:-1755.521   1st Qu.:10128.318   1st Qu.:-1110.456  
 Median :   -0.644   Median :10857.506   Median : -400.623  
 Mean   :    0.000   Mean   :10612.677   Mean   : -172.452  
 3rd Qu.: 1285.666   3rd Qu.:11340.287   3rd Qu.: 1214.203  
 Max.   : 4774.143   Max.   :11502.819   Max.   : 9935.006  
 IQR:
     STL.seasonal STL.trend STL.remainder data 
     3041         1212      2325          3550 
   %  85.7         34.1      65.5         100.0

 Weights:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.0000  0.8385  0.9452  0.8182  0.9845  1.0000 

 Other components: List of 5
 $ win  : Named num [1:3] 931 31 31
 $ deg  : Named int [1:3] 0 1 1
 $ jump : Named num [1:3] 94 4 4
 $ inner: int 1
 $ outer: int 15

La componente estacional tiene mucho peso.

Componentes_SM3_DiariaEneroSTL<- tsSM3_DiariaEnero %>%
  stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_DiariaEneroSTL, ts.colour = 'blue')

Se observa un aumento del consumo energético del lavadero a lo largo de los meses de Enero desde los 8000 vatios hasta los 11000.

plot.ts(tsSM3_DiariaEnero, col = 'gray', main="EvoluciĂ³n del consumo energĂ©tico del AC y termo \n para el mes de Enero")
lines(Componentes_SM3_DiariaEneroSTL$time.series[,2], 
      col = "red", lwd = 1, lty = 2)

PredicciĂ³n con Holt-Winters (suavizado exponencial) para el año 2010

EvoluciĂ³n mensual del consumo (sin estacionalidad). Años 2007-2009

Cocina
# LE quito la estacionalidad
tsSM1_DiariaEnero_detrended <- tsSM1_DiariaEnero - Componentes_SM1_DiariaEneroSTL$time.series[,1]
plot.ts(tsSM1_DiariaEnero_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'EvoluciĂ³n diaria del consumo energĂ©tico de la cocina \n En el mes de Enero',
        cex.main = 0.85)

Vamos a estudiar los resĂ­duos:

par(mfrow = c(1,2))
  plot(Componentes_SM1_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM1_DiariaEneroSTL$time.series[,3])

shapiro.test(Componentes_SM1_DiariaEneroSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM1_DiariaEneroSTL$time.series[, 3]
W = 0.78105, p-value = 1.895e-10

Los resĂ­duos no son solo ruido. AĂºn hay estacionalidad. Volvemos a descomponer la serie desestacionalizada una vez.

## Volvemos a descomponer la serie
plot(stl(tsSM1_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE, t.window = 31))

Parece que la componente irregular no sigue ningĂºn patrĂ³n.

Suavizado exponencial de HoltWinters
  • Primera interacciĂ³n
tsSM1_HW_DiariaEnero<- HoltWinters(tsSM1_DiariaEnero_detrended ,
                                   alpha=TRUE)
summary(tsSM1_HW_DiariaEnero$fitted)
      xhat             level           trend           season        
 Min.   :-8059.6   Min.   :-8138   Min.   :1.099   Min.   :-1783.95  
 1st Qu.:  622.2   1st Qu.: 1337   1st Qu.:1.099   1st Qu.:-1110.06  
 Median : 1780.5   Median : 2239   Median :1.099   Median : -715.70  
 Mean   : 2059.4   Mean   : 2058   Mean   :1.099   Mean   :    0.00  
 3rd Qu.: 3900.9   3rd Qu.: 3276   3rd Qu.:1.099   3rd Qu.:  -41.49  
 Max.   :10897.9   Max.   : 7766   Max.   :1.099   Max.   : 9121.65  
# optimiza
plot(tsSM1_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronĂ³stico de consumo variarĂ¡ en el rango: 0-10897 vatios.hora.

  • Segunda interacciĂ³n:
tsSM1_HW_DiariaEnero<- HoltWinters(tsSM1_DiariaEnero_detrended ,
                                   alpha=TRUE, beta=TRUE)
summary(tsSM1_HW_DiariaEnero$fitted)
      xhat              level           trend               season        
 Min.   :-15283.3   Min.   :-8138   Min.   :-9059.226   Min.   :-1783.95  
 1st Qu.:   247.6   1st Qu.: 1337   1st Qu.:-1401.153   1st Qu.:-1110.06  
 Median :  1412.9   Median : 2239   Median :   -5.467   Median : -715.70  
 Mean   :  2074.9   Mean   : 2058   Mean   :   16.639   Mean   :    0.00  
 3rd Qu.:  4945.5   3rd Qu.: 3276   3rd Qu.: 1526.117   3rd Qu.:  -41.49  
 Max.   : 18063.3   Max.   : 7766   Max.   :13976.323   Max.   : 9121.65  
# optimiza
plot(tsSM1_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Tendencia del consumo creciente.

  • Tercera iteracciĂ³n
tsSM1_HW_DiariaEnero<- HoltWinters(tsSM1_DiariaEnero_detrended ,
                                   alpha=TRUE,
                                   beta=TRUE, gamma=TRUE)
summary(tsSM1_HW_DiariaEnero$fitted)
      xhat              level           trend               season        
 Min.   :-15283.3   Min.   :-8138   Min.   :-9059.226   Min.   :-1783.95  
 1st Qu.:   247.6   1st Qu.: 1337   1st Qu.:-1401.153   1st Qu.:-1110.06  
 Median :  1412.9   Median : 2239   Median :   -5.467   Median : -715.70  
 Mean   :  2074.9   Mean   : 2058   Mean   :   16.639   Mean   :    0.00  
 3rd Qu.:  4945.5   3rd Qu.: 3276   3rd Qu.: 1526.117   3rd Qu.:  -41.49  
 Max.   : 18063.3   Max.   : 7766   Max.   :13976.323   Max.   : 9121.65  
# optimiza
plot(tsSM1_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM1_HW_DiariaEnero_for <- forecast(tsSM1_HW_DiariaEnero, h=31)
plot(tsSM1_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha",
     main="PredicciĂ³n del consumo energĂ©tico de la cocina \n Enero Año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_DiariaEnero_forC <- forecast(tsSM1_HW_DiariaEnero, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha", main="PredicciĂ³n del consumo energĂ©tico de la cocina \n  Enero del 2010", start(2010))

Lavadero
tsSM2_DiariaEnero_detrended <- tsSM2_DiariaEnero - Componentes_SM2_DiariaEneroSTL$time.series[,1]
plot.ts(tsSM2_DiariaEnero_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'EvoluciĂ³n ',
        cex.main = 0.85)

Estudio de los resĂ­duos.

par(mfrow = c(1,2))
  plot(Componentes_SM2_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM2_DiariaEneroSTL$time.series[,3])

shapiro.test(Componentes_SM2_DiariaEneroSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM2_DiariaEneroSTL$time.series[, 3]
W = 0.86799, p-value = 1.342e-07

Los resĂ­duos no siguen una distribuciĂ³n normal. La descomposiciĂ³n no es buena, por lo que volvemos a descomponer la serie.

## Volvemos a descomponer la serie
plot(stl(tsSM2_DiariaEnero_detrended ,s.window = "periodic",
         robust = TRUE,t.window = 31))

Suavizado exponencial de HoltWinters
  • Primera interacciĂ³n
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended ,
                              alpha=TRUE)

summary(tsSM2_HW_DiariaEnero$fitted)
      xhat             level             trend            season       
 Min.   :-5033.6   Min.   :-4350.7   Min.   :-19.53   Min.   :-4239.4  
 1st Qu.:  109.6   1st Qu.:  551.7   1st Qu.:-19.53   1st Qu.:-1085.6  
 Median : 1615.7   Median : 2193.8   Median :-19.53   Median : -253.8  
 Mean   : 2167.6   Mean   : 2187.1   Mean   :-19.53   Mean   :    0.0  
 3rd Qu.: 4042.5   3rd Qu.: 2820.3   3rd Qu.:-19.53   3rd Qu.:  903.3  
 Max.   :11795.2   Max.   : 9030.9   Max.   :-19.53   Max.   : 6231.5  
plot(tsSM2_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Rango de pronĂ³stico del consumo energĂ©tico del lavadero: 0-11795 vatios-hora.

  • Segunda interacciĂ³n: tendencia
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended ,
                              alpha=TRUE,beta=TRUE)

summary(tsSM2_HW_DiariaEnero$fitted)
      xhat            level             trend                season       
 Min.   :-11600   Min.   :-4350.7   Min.   :-10691.968   Min.   :-4239.4  
 1st Qu.: -2922   1st Qu.:  551.7   1st Qu.: -2362.484   1st Qu.:-1085.6  
 Median :  1808   Median : 2193.8   Median :    -8.537   Median : -253.8  
 Mean   :  2159   Mean   : 2187.1   Mean   :   -27.645   Mean   :    0.0  
 3rd Qu.:  6071   3rd Qu.: 2820.3   3rd Qu.:  2539.702   3rd Qu.:  903.3  
 Max.   : 19223   Max.   : 9030.9   Max.   : 10338.871   Max.   : 6231.5  
plot(tsSM2_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

Aumento del consumo

  • Tercera interacciĂ³n
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended ,
                              alpha=TRUE,
                              beta=TRUE, gamma=TRUE)

summary(tsSM2_HW_DiariaEnero$fitted)
      xhat            level             trend                season       
 Min.   :-11600   Min.   :-4350.7   Min.   :-10691.968   Min.   :-4239.4  
 1st Qu.: -2922   1st Qu.:  551.7   1st Qu.: -2362.484   1st Qu.:-1085.6  
 Median :  1808   Median : 2193.8   Median :    -8.537   Median : -253.8  
 Mean   :  2159   Mean   : 2187.1   Mean   :   -27.645   Mean   :    0.0  
 3rd Qu.:  6071   3rd Qu.: 2820.3   3rd Qu.:  2539.702   3rd Qu.:  903.3  
 Max.   : 19223   Max.   : 9030.9   Max.   : 10338.871   Max.   : 6231.5  
plot(tsSM2_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones, es decir, el consumo energĂ©tico del mes de enero del año 2010.

## HoltWinters forecast & plot
tsSM2_HW_DiariaEnero_for <- forecast(tsSM2_HW_DiariaEnero, h=31)
plot(tsSM2_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n del consumo energĂ©tico del lavadero \n Enero del año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_DiariaEnero_forC <- forecast(tsSM2_HW_DiariaEnero, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n del consumo energĂ©tico del lavadero \n Enero del año 2010", start(2010))

AC y termo eléctrico
tsSM3_DiariaEnero_detrended <- tsSM3_DiariaEnero - Componentes_SM3_DiariaEneroSTL$time.series[,1]
plot.ts(tsSM3_DiariaEnero_detrended, 
        ylab = "EnergĂ­a (Vatios-hora)", 
        main = 'VariaciĂ³n el consumo energĂ©tico del mes de Enero',
        cex.main = 0.85)

Vamos a estudiar los resĂ­duos:

par(mfrow = c(1,2))
  plot(Componentes_SM3_DiariaEneroSTL$time.series[,3], col = "blue",
       main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM3_DiariaEneroSTL$time.series[,3])

shapiro.test(Componentes_SM3_DiariaEneroSTL$time.series[,3])

    Shapiro-Wilk normality test

data:  Componentes_SM3_DiariaEneroSTL$time.series[, 3]
W = 0.90263, p-value = 3.877e-06

Los resĂ­duos no siguen una distribuciĂ³n normal. La descomposiciĂ³n no es buena, por lo que vamos a volver a descomponer la serie:

## Volvemos a descomponer la serie
plot(stl(tsSM3_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE,t.window = 31))

La componente irregular no presenta ningĂºn patrĂ³n. AdemĂ¡s, observamos una fuerte tendencia creciente del consumo.

Suavizado exponencial de HoltWinters
  • Primera interacciĂ³n: rango del pronĂ³stico de consumo.
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_DiariaEnero <- HoltWinters(tsSM3_DiariaEnero_detrended ,
                             alpha=TRUE)
summary(tsSM3_HW_DiariaEnero$fitted)
      xhat           level             trend            season       
 Min.   :-1411   Min.   : -515.4   Min.   :-29.16   Min.   :-7769.7  
 1st Qu.: 8391   1st Qu.: 8869.1   1st Qu.:-29.16   1st Qu.: -971.9  
 Median :10175   Median :10509.4   Median :-29.16   Median : -245.3  
 Mean   :10326   Mean   :10354.9   Mean   :-29.16   Mean   :    0.0  
 3rd Qu.:12873   3rd Qu.:12369.5   3rd Qu.:-29.16   3rd Qu.:  810.6  
 Max.   :22353   Max.   :19582.7   Max.   :-29.16   Max.   : 4510.7  
plot(tsSM3_HW_DiariaEnero,
     xlab="Fecha",
     ylab="Valores observados y ajustados",
     main="Suavizado exponencial de Holt-Winters")

El pronĂ³stico de consumo energĂ©tico del termo elĂ©ctrico y el aire acondicionado variarĂ¡ en torno al siguiente rango: 0-22353 vatios-hora

PronĂ³stico de HoltWinters para el año 2010

Vamos a predecir las prĂ³ximas 12 observaciones.

## HoltWinters forecast & plot
tsSM3_HW_DiariaEnero_for <- forecast(tsSM3_HW_DiariaEnero, h=31)
plot(tsSM3_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
     main="PredicciĂ³n del consumo energĂ©tico del AC y termo elĂ©ctrico \n Mes de enero del año 2010" ) 

## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_DiariaEnero_forC <- forecast( object = tsSM3_HW_DiariaEnero_for, h=31)
## Plot only the forecasted area
plot(tsSM3_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo elĂ©ctrico", main="PredicciĂ³n para el mes de Enero del año 2010", start(2010), ylim=c(0,15000))

ComparaciĂ³n de resultados. PredicciĂ³n despues de eliminar estacionalidad

seriesDiariasEnero <- Granularidad_dias %>% 
  filter( year == 2010 & month == 1)


RMSE_SM1_Enero<-accuracy(tsSM1_HW_DiariaEnero_forC ,seriesDiariasEnero $Sub_metering_1)[4]
RMSE_SM2_Enero<-accuracy(tsSM2_HW_DiariaEnero_forC  ,seriesDiariasEnero $Sub_metering_2)[4]
RMSE_SM3_Enero<-accuracy(tsSM3_HW_DiariaEnero_forC  ,seriesDiariasEnero $Sub_metering_3)[4]



Error<-
cbind.data.frame(
      Submedidor=c("Cocina","Lavadero","AC y TermoE"),
      RMSE = c(RMSE_SM1_Enero,
               RMSE_SM2_Enero,
               RMSE_SM3_Enero)
   
    
)

Error %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")
Submedidor RMSE
Cocina 150276.386
Lavadero 36593.987
AC y TermoE 7385.421

El error para la cocina es extremadamente alto, sin embargo para el lavadero podría ser aceptable y para el termo eléctrico y aire acondicionado es bastante bajo, el resultado es bueno.

PronĂ³stico del consumo para los siguientes meses de Enero

Cocina

# Creamos la serie temporal con los datos completos
seriesDiariasEnero <- Granularidad_dias %>% filter( month == 1)
## Creamos un objeto TS
   tsSM1_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_1, frequency=31, start=c(2007,1))
   
   tsSM2_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_2, frequency=31, start=c(2007,1))
   
   tsSM3_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_3, frequency=31, start=c(2007,1))
# Descomponemos la serie
Componentes_SM1_DiariaEneroSTL<-
  tsSM1_DiariaEnero %>%stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM1_DiariaEneroSTL, ts.colour = 'blue')

# EliminaciĂ³n de la estacionalidad
tsSM1_DiariaEnero_detrended <- tsSM1_DiariaEnero - Componentes_SM1_DiariaEneroSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM1_HW_DiariaEnero <- HoltWinters(tsSM1_DiariaEnero_detrended   , 
                                seasonal = "additive",
  beta=TRUE, 
  gamma=TRUE,
   alpha=TRUE
)
# PronĂ³stico para los siguientes 13 meses 
PrediccionFinalCocinaEnero <- forecast(tsSM1_HW_DiariaEnero, h=31)
## Plot only the forecasted area
plot(PrediccionFinalCocinaEnero, ylab= "Vatios-Hora", xlab="Fecha", main="PredicciĂ³n el consumo energĂ©tico \n  de la cocina para Enero del 2011", start(2011), ylim = c(0,5000))

Las predicciones son las siguientes:

PrediccionFinalCocinaEnero[["mean"]] 
Time Series:
Start = c(2011, 1) 
End = c(2011, 31) 
Frequency = 31 
 [1]   3817.7251   1703.2496   1798.5550   1789.3679   1740.9831   1714.8134
 [7]    191.6482   1149.2603   -742.6860  -1190.2384  -1802.0016   4258.6795
[13]  -2247.3347    793.1571  -3181.4382  -3742.3316  -2512.3870  -3889.7260
[19]  -4244.3764   -788.0276   4649.3493  -4534.0677  -7381.3916  -5918.3862
[25]  -6165.6491  -7019.4282  -6945.7927  -6447.0977  -8353.3125  -9155.1739
[31] -10117.8606

Vimos que estas predicciones tenĂ­an un error muy grande por lo que no son muy fiables.

Lavadero

# Descomponemos la serie
Componentes_SM2_DiariaEneroSTL<-
  tsSM2_DiariaEnero %>%stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM2_DiariaEneroSTL, ts.colour = 'blue')

# EliminaciĂ³n de la estacionalidad
tsSM2_DiariaEnero_detrended <- tsSM2_DiariaEnero - Componentes_SM2_DiariaEneroSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended   , 
                                seasonal = "additive",
   beta=TRUE,
 gamma=TRUE
#alpha=TRUE
)
# PronĂ³stico para los siguientes 13 meses 
PrediccionFinalLavaderoEnero <- forecast(tsSM2_HW_DiariaEnero , h=31)
## Plot only the forecasted area
plot(PrediccionFinalLavaderoEnero, ylab= "Vatios-Hora", xlab="Fecha", main="PredicciĂ³n el consumo energĂ©tico \n  del lavadero para Enero del 2011", start(2011), ylim=c(0,10000))

Las predicciones son las siguientes:

PrediccionFinalLavaderoEnero[["mean"]] 
Time Series:
Start = c(2011, 1) 
End = c(2011, 31) 
Frequency = 31 
 [1] -3208.6987 -2451.0215   337.0196 -8998.2347 -3419.9241 -2154.1838
 [7] -4699.1157 -4641.6730 -2098.9702 -5915.8668 -5908.5963 -3312.5989
[13] -4924.5015 -5895.6914 -5117.6370 -7699.2853 -9518.6745 -7476.8637
[19] -6785.2376 -6129.8531 -8420.1070 -7299.7906 -3701.8091 -6040.7288
[25] -7146.2346 -5477.4638 -2360.3339 -6951.1231 -4149.3855  2215.0182
[31] -7088.5707

Predicciones negativas.

Aire acondicionado y termo eléctrico

# Descomponemos la serie
Componentes_SM3_DiariaEneroSTL<-
  tsSM3_DiariaEnero %>%stl( t.window=31,
  s.window="periodic",  robust=TRUE)
autoplot(Componentes_SM3_DiariaEneroSTL, ts.colour = 'blue')

# EliminaciĂ³n de la estacionalidad
tsSM3_DiariaEnero_detrended <- tsSM3_DiariaEnero - Componentes_SM3_DiariaEneroSTL$time.series[,1]

# Suavizado exponencial tal y como hemos obtenido el mejor modelo
tsSM3_HW_DiariaEnero <- HoltWinters(tsSM3_DiariaEnero_detrended   , 
        seasonal = "additive",
 #  beta=TRUE
  #gamma=TRUE
  alpha=TRUE
)
plot(tsSM3_HW_DiariaEnero)

# PronĂ³stico para los siguientes 13 meses 
PrediccionFinalACTermoEnero <- forecast(tsSM3_HW_DiariaEnero , h=31)
## Plot only the forecasted area
plot(PrediccionFinalACTermoEnero, ylab= "Vatios-Hora", xlab="Fecha", main="PredicciĂ³n el consumo energĂ©tico \n  del AC y Termo para Enero del 2011", start(2011), ylim = c(0,20000))

Las predicciones son las siguientes:

PrediccionFinalACTermoEnero[["mean"]] 
Time Series:
Start = c(2011, 1) 
End = c(2011, 31) 
Frequency = 31 
 [1] 17286.890  9397.946  7612.869 10384.247 10971.474  8229.004  7210.607
 [8]  8490.469  7574.731  6991.382  8002.352  7186.606  6843.313  6456.335
[15]  9699.005  6319.642 10407.530  9105.477  4273.980  9062.211  7955.532
[22]  9057.087  1206.155 10639.813  7764.096  5200.978 10242.690  7921.854
[29]  9748.997  6640.078 11432.974

Predicciones negativas.